Code and Simulations for: Plastic responses to novel environments are biased towards dimensions with high genetic variance

1.0 – Simulation 1: G and P Comparisons Within and Across Environments

1.1 – Setting up and Loading the Packages and Data

1.1.1 – Packages and Functions

Lets first load the relevant packages and functions. Note that, because the .Rmd file is located within the R/ folder we need to modify the source script for the “homemade” functions. Note that we have # out kableExtra as this only needs to be done once.

    rm(list=ls())
    source("./R/GP_compare_func.R")
    source("./R/evolvability_func.R")
    pacman::p_load(grid, png, shape, car, tidyr, dplyr, kableExtra, latex2exp, metafor, brms, rstan, timetree, MCMCglmm, data.table, ggplot2, knitr, rmarkdown)
    #devtools::install_github("haozhu233/kableExtra")

1.1.2 – Loading the Matrices from the data/ Folder

We want to be able to get into all the paper folders, bring in all the matrices and information about traits, perform operations on the matrices and extract the relevant information from them so we can do Monte Carlo simulations. To start, we can acquire the path names of all the matrix and description files within each paper folder of the matrices in data_GP_compare/.The below code will use a function that is custom built to grab all the above information.

      matrices <- files(dir = "./data_GP_compare/", type = "matrix", file_pattern = "[WRs][BR]")
    descripDat <- files(dir = "./data_GP_compare/", type = "D", file_pattern = "[WRs][BR]")

1.1.3 – Load the Meta-Data

The meta data is needed as it contains the relevant moderators and information for each of the papers that are included in order to merge this data with the simulations before running the meta-regression models. We can also see the data in the table below:

    meta <- read.csv("./data_GP_compare/stdy_metadata_final.csv", stringsAsFactor = FALSE)
    kable(meta) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "200px")

file	authors	year	study	env	type	n	avg_n_env	avg_n_type	trait_num	col	genus	species	design_raw	design	group	name.type	nam.env	env.comp
./data_GP_compare/RR073/B_n=68_G.csv	Evans et al.	2015	RR073	B	G	68	114.5	66.0	8	2	Poecilia	reticulata	half-sib	half-sib	fish	RR073.G	RR073.B	BS
./data_GP_compare/RR073/B_n=161_P.csv	Evans et al.	2015	RR073	B	P	161	114.5	153.5	8	1	Poecilia	reticulata	half-sib	half-sib	fish	RR073.P	RR073.B	BS
./data_GP_compare/RR073/S1_n=64_G.csv	Evans et al.	2015	RR073	S	G	64	105.0	66.0	8	4	Poecilia	reticulata	half-sib	half-sib	fish	RR073.G	RR073.S	BS
./data_GP_compare/RR073/S1_n=146_P.csv	Evans et al.	2015	RR073	S	P	146	105.0	153.5	8	3	Poecilia	reticulata	half-sib	half-sib	fish	RR073.P	RR073.S	BS
./data_GP_compare/RR076/B_n=69_G.csv	Guan et al.	2016	RR076	B	G	69	1097.0	69.0	3	6	Scophthalmus	maximus	half-sib	half-sib	fish	RR076.G	RR076.B	BS
./data_GP_compare/RR076/B_n=2125_P.csv	Guan et al.	2016	RR076	B	P	2125	1097.0	2525.0	3	5	Scophthalmus	maximus	half-sib	half-sib	fish	RR076.P	RR076.B	BS
./data_GP_compare/RR076/S1_n=69_G.csv	Guan et al.	2016	RR076	S	G	69	1497.0	69.0	3	8	Scophthalmus	maximus	half-sib	half-sib	fish	RR076.G	RR076.S	BS
./data_GP_compare/RR076/S1_n=2925_P.csv	Guan et al.	2016	RR076	S	P	2925	1497.0	2525.0	3	7	Scophthalmus	maximus	half-sib	half-sib	fish	RR076.P	RR076.S	BS
./data_GP_compare/RR087a/A_n=100_G.csv	Sae-Lim et al.	2013	RR087a	A	G	100	1232.0	100.0	3	9	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087a.G	RR087a.A	AN
./data_GP_compare/RR087a/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087a	A	P	2364	1232.0	2127.5	3	10	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087a.P	RR087a.A	AN
./data_GP_compare/RR087a/N1_n=100_G.csv	Sae-Lim et al.	2013	RR087a	N	G	100	995.5	100.0	3	11	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087a.G	RR087a.N	AN
./data_GP_compare/RR087a/N1_n=1891_P.csv	Sae-Lim et al.	2013	RR087a	N	P	1891	995.5	2127.5	3	12	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087a.P	RR087a.N	AN
./data_GP_compare/RR087b/A_n=100_G.csv	Sae-Lim et al.	2013	RR087b	A	G	100	1232.0	100.0	3	13	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087b.G	RR087b.A	AN
./data_GP_compare/RR087b/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087b	A	P	2364	1232.0	2577.0	3	14	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087b.P	RR087b.A	AN
./data_GP_compare/RR087b/N2_n=100_G.csv	Sae-Lim et al.	2013	RR087b	N	G	100	1445.0	100.0	3	15	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087b.G	RR087b.N	AN
./data_GP_compare/RR087b/N2_n=2790_P.csv	Sae-Lim et al.	2013	RR087b	N	P	2790	1445.0	2577.0	3	16	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087b.P	RR087b.N	AN
./data_GP_compare/RR087c/A_n=100_G.csv	Sae-Lim et al.	2013	RR087c	A	G	100	1232.0	100.0	3	17	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087c.G	RR087c.A	AN
./data_GP_compare/RR087c/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087c	A	P	2364	1232.0	2091.0	3	18	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087c.P	RR087c.A	AN
./data_GP_compare/RR087c/N3_n=100_G.csv	Sae-Lim et al.	2013	RR087c	N	G	100	959.0	100.0	3	19	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087c.G	RR087c.N	AN
./data_GP_compare/RR087c/N3_n=1818_P.csv	Sae-Lim et al.	2013	RR087c	N	P	1818	959.0	2091.0	3	20	Oncorhynchus	mykiss	full-sib	full-sib	fish	RR087c.P	RR087c.N	AN
./data_GP_compare/sR097/B_n=26_G.csv	Kasule	1991	sR097	B	G	26	169.0	26.0	2	21	Dysdercus	fasciatus	half-sib	half-sib	insect	sR097.G	sR097.B	BS
./data_GP_compare/sR097/B_n=312_P.csv	Kasule	1991	sR097	B	P	312	169.0	312.0	2	22	Dysdercus	fasciatus	half-sib	half-sib	insect	sR097.P	sR097.B	BS
./data_GP_compare/sR097/S1_n=26_G.csv	Kasule	1991	sR097	S	G	26	169.0	26.0	2	23	Dysdercus	fasciatus	half-sib	half-sib	insect	sR097.G	sR097.S	BS
./data_GP_compare/sR097/S1_n=312_P.csv	Kasule	1991	sR097	S	P	312	169.0	312.0	2	24	Dysdercus	fasciatus	half-sib	half-sib	insect	sR097.P	sR097.S	BS
./data_GP_compare/sR110b/B_n=37_G.csv	Blanckenhorn & Heyland	2004	sR110b	B	G	37	82.5	32.0	2	26	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110b.G	sR110b.B	BS
./data_GP_compare/sR110b/B_n=128_P.csv	Blanckenhorn & Heyland	2004	sR110b	B	P	128	82.5	113.0	2	25	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110b.P	sR110b.B	BS
./data_GP_compare/sR110b/S1_n=27_G.csv	Blanckenhorn & Heyland	2004	sR110b	S	G	27	62.5	32.0	2	27	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110b.G	sR110b.S	BS
./data_GP_compare/sR110b/S1_n=98_P.csv	Blanckenhorn & Heyland	2004	sR110b	S	P	98	62.5	113.0	2	28	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110b.P	sR110b.S	BS
./data_GP_compare/sR110c/B_n=40_G.csv	Blanckenhorn & Heyland	2004	sR110c	B	G	40	105.5	36.5	2	30	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110c.G	sR110c.B	BS
./data_GP_compare/sR110c/B_n=171_P.csv	Blanckenhorn & Heyland	2004	sR110c	B	P	171	105.5	150.0	2	29	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110c.P	sR110c.B	BS
./data_GP_compare/sR110c/S1_n=33_G.csv	Blanckenhorn & Heyland	2004	sR110c	S	G	33	81.0	36.5	2	32	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110c.G	sR110c.S	BS
./data_GP_compare/sR110c/S1_n=129_P.csv	Blanckenhorn & Heyland	2004	sR110c	S	P	129	81.0	150.0	2	31	Scathophaga	stercoraria	half-sib	half-sib	insect	sR110c.P	sR110c.S	BS
./data_GP_compare/sR115/A_n=350_G.csv	Conner et al.	2003	sR115	A	G	350	741.5	471.5	6	34	Raphanus	raphanistrum	half-sib	half-sib	plant	sR115.G	sR115.A	AN
./data_GP_compare/sR115/A_n=1133_P.csv	Conner et al.	2003	sR115	A	P	1133	741.5	1010.5	6	33	Raphanus	raphanistrum	half-sib	half-sib	plant	sR115.P	sR115.A	AN
./data_GP_compare/sR115/N_n=593_G.csv	Conner et al.	2003	sR115	N	G	593	740.5	471.5	6	35	Raphanus	raphanistrum	half-sib	half-sib	plant	sR115.G	sR115.N	AN
./data_GP_compare/sR115/N_n=888_P.csv	Conner et al.	2003	sR115	N	P	888	740.5	1010.5	6	36	Raphanus	raphanistrum	half-sib	half-sib	plant	sR115.P	sR115.N	AN
./data_GP_compare/sR116/A_n=366_G.csv	Delcourt et al.	2010	sR116	A	G	366	1840.5	365.5	8	38	Drosophila	serrata	half-sib	half-sib	insect	sR116.G	sR116.A	AN
./data_GP_compare/sR116/A_n=3315_P.csv	Delcourt et al.	2010	sR116	A	P	3315	1840.5	3299.0	8	37	Drosophila	serrata	half-sib	half-sib	insect	sR116.P	sR116.A	AN
./data_GP_compare/sR116/N1_n=365_G.csv	Delcourt et al.	2010	sR116	N	G	365	1824.0	365.5	8	40	Drosophila	serrata	half-sib	half-sib	insect	sR116.G	sR116.N	AN
./data_GP_compare/sR116/N1_n=3283_P.csv	Delcourt et al.	2010	sR116	N	P	3283	1824.0	3299.0	8	39	Drosophila	serrata	half-sib	half-sib	insect	sR116.P	sR116.N	AN
./data_GP_compare/sR145/B_n=52_G.csv	Tucic et al.	1991	sR145	B	G	52	101.0	47.5	8	42	Acanthoscelides	obtectus	full-sib	full-sib	insect	sR145.G	sR145.B	BS
./data_GP_compare/sR145/B_n=150_P.csv	Tucic et al.	1991	sR145	B	P	150	101.0	137.0	8	41	Acanthoscelides	obtectus	full-sib	full-sib	insect	sR145.P	sR145.B	BS
./data_GP_compare/sR145/S1_n=43_G.csv	Tucic et al.	1991	sR145	S	G	43	83.5	47.5	8	44	Acanthoscelides	obtectus	full-sib	full-sib	insect	sR145.G	sR145.S	BS
./data_GP_compare/sR145/S1_n=124_P.csv	Tucic et al.	1991	sR145	S	P	124	83.5	137.0	8	43	Acanthoscelides	obtectus	full-sib	full-sib	insect	sR145.P	sR145.S	BS
./data_GP_compare/WB002a/A_n=39_G.csv	Begin & Roff	2001	WB002a	A	G	39	372.0	39.0	5	45	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB002a.G	WB002a.A	AN
./data_GP_compare/WB002a/A_n=705_P.csv	Begin & Roff	2001	WB002a	A	P	705	372.0	605.0	5	46	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB002a.P	WB002a.A	AN
./data_GP_compare/WB002a/N1_n=39_G.csv	Begin & Roff	2001	WB002a	N	G	39	272.0	39.0	5	47	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB002a.G	WB002a.N	AN
./data_GP_compare/WB002a/N1_n=505_P.csv	Begin & Roff	2001	WB002a	N	P	505	272.0	605.0	5	48	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB002a.P	WB002a.N	AN
./data_GP_compare/WB003/A_n=56_G.csv	Begin et al.	2004	WB003	A	G	56	294.5	59.0	5	50	Gryllus	firmus	full-sib	full-sib	insect	WB003.G	WB003.A	AN
./data_GP_compare/WB003/A_n=533_P.csv	Begin et al.	2004	WB003	A	P	533	294.5	697.5	5	49	Gryllus	firmus	full-sib	full-sib	insect	WB003.P	WB003.A	AN
./data_GP_compare/WB003/N1_n=62_G.csv	Begin et al.	2004	WB003	N	G	62	387.5	59.0	5	51	Gryllus	firmus	full-sib	full-sib	insect	WB003.G	WB003.N	AN
./data_GP_compare/WB003/N1_n=862_P.csv	Begin et al.	2004	WB003	N	P	862	387.5	697.5	5	52	Gryllus	firmus	full-sib	full-sib	insect	WB003.P	WB003.N	AN
./data_GP_compare/WB003_2/A_n=56_G.csv	Begin et al.	2004	WB003_2	A	G	56	294.5	58.0	5	54	Gryllus	firmus	full-sib	full-sib	insect	WB003_2.G	WB003_2.A	AN
./data_GP_compare/WB003_2/A_n=533_P.csv	Begin et al.	2004	WB003_2	A	P	533	294.5	624.0	5	53	Gryllus	firmus	full-sib	full-sib	insect	WB003_2.P	WB003_2.A	AN
./data_GP_compare/WB003_2/N2_n=60_G.csv	Begin et al.	2004	WB003_2	N	G	60	462.0	58.0	5	55	Gryllus	firmus	full-sib	full-sib	insect	WB003_2.G	WB003_2.N	AN
./data_GP_compare/WB003_2/N2_n=715_P.csv	Begin et al.	2004	WB003_2	N	P	715	462.0	624.0	5	56	Gryllus	firmus	full-sib	full-sib	insect	WB003_2.P	WB003_2.N	AN
./data_GP_compare/WB003a/A_n=45_G.csv	Begin et al.	2004	WB003a	A	G	45	153.5	43.0	5	58	Gryllus	firmus	full-sib	full-sib	insect	WB003a.G	WB003a.A	AN
./data_GP_compare/WB003a/A_n=262_P.csv	Begin et al.	2004	WB003a	A	P	262	153.5	288.5	5	57	Gryllus	firmus	full-sib	full-sib	insect	WB003a.P	WB003a.A	AN
./data_GP_compare/WB003a/N1_n=41_G.csv	Begin et al.	2004	WB003a	N	G	41	82.5	43.0	5	60	Gryllus	firmus	full-sib	full-sib	insect	WB003a.G	WB003a.N	AN
./data_GP_compare/WB003a/N1_n=315_P.csv	Begin et al.	2004	WB003a	N	P	315	82.5	288.5	5	59	Gryllus	firmus	full-sib	full-sib	insect	WB003a.P	WB003a.N	AN
./data_GP_compare/WB003a_2/A_n=45_G.csv	Begin et al.	2004	WB003a_2	A	G	45	153.5	38.5	5	62	Gryllus	firmus	full-sib	full-sib	insect	WB003a_2.G	WB003a_2.A	AN
./data_GP_compare/WB003a_2/A_n=262_P.csv	Begin et al.	2004	WB003a_2	A	P	262	153.5	197.5	5	61	Gryllus	firmus	full-sib	full-sib	insect	WB003a_2.P	WB003a_2.A	AN
./data_GP_compare/WB003a_2/N2_n=32_G.csv	Begin et al.	2004	WB003a_2	N	G	32	178.0	38.5	5	64	Gryllus	firmus	full-sib	full-sib	insect	WB003a_2.G	WB003a_2.N	AN
./data_GP_compare/WB003a_2/N2_n=133_P.csv	Begin et al.	2004	WB003a_2	N	P	133	178.0	197.5	5	63	Gryllus	firmus	full-sib	full-sib	insect	WB003a_2.P	WB003a_2.N	AN
./data_GP_compare/WB004/A_n=300_G.csv	Bennington & McGraw	2000	WB004	A	G	300	342.0	300.0	5	65	Impatiens	pallida	half-sib	half-sib	plant	WB004.G	WB004.A	AN
./data_GP_compare/WB004/A_n=384_P.csv	Bennington & McGraw	2000	WB004	A	P	384	342.0	244.5	5	66	Impatiens	pallida	half-sib	half-sib	plant	WB004.P	WB004.A	AN
./data_GP_compare/WB004/N1_n=300_G.csv	Bennington & McGraw	2000	WB004	N	G	300	202.5	300.0	5	68	Impatiens	pallida	half-sib	half-sib	plant	WB004.G	WB004.N	AN
./data_GP_compare/WB004/N1_n=105_P.csv	Bennington & McGraw	2000	WB004	N	P	105	202.5	244.5	5	67	Impatiens	pallida	half-sib	half-sib	plant	WB004.P	WB004.N	AN
./data_GP_compare/WB004a/A_n=300_G.csv	Bennington & McGraw	2000	WB004a	A	G	300	190.0	300.0	5	69	Impatiens	pallida	half-sib	half-sib	plant	WB004a.G	WB004a.A	AN
./data_GP_compare/WB004a/A_n=80_P.csv	Bennington & McGraw	2000	WB004a	A	P	80	190.0	108.5	5	70	Impatiens	pallida	half-sib	half-sib	plant	WB004a.P	WB004a.A	AN
./data_GP_compare/WB004a/N1_n=300_G.csv	Bennington & McGraw	2000	WB004a	N	G	300	218.5	300.0	5	72	Impatiens	pallida	half-sib	half-sib	plant	WB004a.G	WB004a.N	AN
./data_GP_compare/WB004a/N1_n=137_P.csv	Bennington & McGraw	2000	WB004a	N	P	137	218.5	108.5	5	71	Impatiens	pallida	half-sib	half-sib	plant	WB004a.P	WB004a.N	AN
./data_GP_compare/WB007/B_n=132_G.csv	Brock et al.	2010	WB007	B	G	132	366.5	131.5	11	73	Brassica	rapa	RIL	half-sib	plant	WB007.G	WB007.B	BS
./data_GP_compare/WB007/B_n=601_P.csv	Brock et al.	2010	WB007	B	P	601	366.5	562.5	11	74	Brassica	rapa	RIL	half-sib	plant	WB007.P	WB007.B	BS
./data_GP_compare/WB007/S1_n=131_G.csv	Brock et al.	2010	WB007	S	G	131	327.5	131.5	11	75	Brassica	rapa	RIL	half-sib	plant	WB007.G	WB007.S	BS
./data_GP_compare/WB007/S1_n=524_P.csv	Brock et al.	2010	WB007	S	P	524	327.5	562.5	11	76	Brassica	rapa	RIL	half-sib	plant	WB007.P	WB007.S	BS
./data_GP_compare/WB010/A_n=141_G.csv	Collins et al.	1999	WB010	A	G	141	840.5	101.5	4	77	Achroia	grisella	half-sib	half-sib	insect	WB010.G	WB010.A	AN
./data_GP_compare/WB010/A_n=1540_P.csv	Collins et al.	1999	WB010	A	P	1540	840.5	991.5	4	78	Achroia	grisella	half-sib	half-sib	insect	WB010.P	WB010.A	AN
./data_GP_compare/WB010/N1_n=62_G.csv	Collins et al.	1999	WB010	N	G	62	252.5	101.5	4	80	Achroia	grisella	half-sib	half-sib	insect	WB010.G	WB010.N	AN
./data_GP_compare/WB010/N1_n=443_P.csv	Collins et al.	1999	WB010	N	P	443	252.5	991.5	4	79	Achroia	grisella	half-sib	half-sib	insect	WB010.P	WB010.N	AN
./data_GP_compare/WB011/B_n=404_G.csv	Czesak & Fox	2003	WB011	B	G	404	1535.5	404.0	2	82	Stator	limbatus	half-sib	half-sib	insect	WB011.G	WB011.B	BS
./data_GP_compare/WB011/B_n=2667_P.csv	Czesak & Fox	2003	WB011	B	P	2667	1535.5	2667.0	2	81	Stator	limbatus	half-sib	half-sib	insect	WB011.P	WB011.B	BS
./data_GP_compare/WB011/S1_n=404_G.csv	Czesak & Fox	2003	WB011	S	G	404	1535.5	404.0	2	84	Stator	limbatus	half-sib	half-sib	insect	WB011.G	WB011.S	BS
./data_GP_compare/WB011/S1_n=2667_P.csv	Czesak & Fox	2003	WB011	S	P	2667	1535.5	2667.0	2	83	Stator	limbatus	half-sib	half-sib	insect	WB011.P	WB011.S	BS
./data_GP_compare/WB028/B_n=60_G.csv	Klause & Morin	2001	WB028	B	G	60	330.0	60.0	2	85	Priophorus	pallipes	full-sib	full-sib	insect	WB028.G	WB028.B	BS
./data_GP_compare/WB028/B_n=600_P.csv	Klause & Morin	2001	WB028	B	P	600	330.0	600.0	2	86	Priophorus	pallipes	full-sib	full-sib	insect	WB028.P	WB028.B	BS
./data_GP_compare/WB028/S1_n=60_G.csv	Klause & Morin	2001	WB028	S	G	60	330.0	60.0	2	87	Priophorus	pallipes	full-sib	full-sib	insect	WB028.G	WB028.S	BS
./data_GP_compare/WB028/S1_n=600_P.csv	Klause & Morin	2001	WB028	S	P	600	330.0	600.0	2	88	Priophorus	pallipes	full-sib	full-sib	insect	WB028.P	WB028.S	BS
./data_GP_compare/WB029/B_n=63_G.csv	King et al.	2011	WB029	B	G	63	862.5	63.0	3	90	Gryllus	firmus	half-sib	half-sib	insect	WB029.G	WB029.B	BS
./data_GP_compare/WB029/B_n=1662_P.csv	King et al.	2011	WB029	B	P	1662	862.5	1695.0	3	89	Gryllus	firmus	half-sib	half-sib	insect	WB029.P	WB029.B	BS
./data_GP_compare/WB029/S1_n=63_G.csv	King et al.	2011	WB029	S	G	63	895.5	63.0	3	92	Gryllus	firmus	half-sib	half-sib	insect	WB029.G	WB029.S	BS
./data_GP_compare/WB029/S1_n=1728_P.csv	King et al.	2011	WB029	S	P	1728	895.5	1695.0	3	91	Gryllus	firmus	half-sib	half-sib	insect	WB029.P	WB029.S	BS
./data_GP_compare/WB033/A_n=60_G.csv	Lau et al.	2014	WB033	A	G	60	293.0	60.0	4	94	Arabidopsis	thaliana	RIL	half-sib	plant	WB033.G	WB033.A	AN
./data_GP_compare/WB033/A_n=526_P.csv	Lau et al.	2014	WB033	A	P	526	293.0	513.0	4	93	Arabidopsis	thaliana	RIL	half-sib	plant	WB033.P	WB033.A	AN
./data_GP_compare/WB033/N1_n=60_G.csv	Lau et al.	2014	WB033	N	G	60	280.0	60.0	4	96	Arabidopsis	thaliana	RIL	half-sib	plant	WB033.G	WB033.N	AN
./data_GP_compare/WB033/N1_n=500_P.csv	Lau et al.	2014	WB033	N	P	500	280.0	513.0	4	95	Arabidopsis	thaliana	RIL	half-sib	plant	WB033.P	WB033.N	AN
./data_GP_compare/WB033a/A_n=60_G.csv	Lau et al.	2014	WB033a	A	G	60	259.5	60.0	4	98	Arabidopsis	thaliana	RIL	half-sib	plant	WB033a.G	WB033a.A	AN
./data_GP_compare/WB033a/A_n=459_P.csv	Lau et al.	2014	WB033a	A	P	459	259.5	464.5	4	97	Arabidopsis	thaliana	RIL	half-sib	plant	WB033a.P	WB033a.A	AN
./data_GP_compare/WB033a/N1_n=60_G.csv	Lau et al.	2014	WB033a	N	G	60	265.0	60.0	4	100	Arabidopsis	thaliana	RIL	half-sib	plant	WB033a.G	WB033a.N	AN
./data_GP_compare/WB033a/N1_n=470_P.csv	Lau et al.	2014	WB033a	N	P	470	265.0	464.5	4	99	Arabidopsis	thaliana	RIL	half-sib	plant	WB033a.P	WB033a.N	AN
./data_GP_compare/WB035/B_n=432_G.csv	Messina et al.	2013	WB035	B	G	432	1316.5	432.0	3	102	Callosobruchus	maculatus	half-sib	half-sib	insect	WB035.G	WB035.B	BS
./data_GP_compare/WB035/B_n=2201_P.csv	Messina et al.	2013	WB035	B	P	2201	1316.5	2204.0	3	101	Callosobruchus	maculatus	half-sib	half-sib	insect	WB035.P	WB035.B	BS
./data_GP_compare/WB035/S1_n=432_G.csv	Messina et al.	2013	WB035	S	G	432	1319.5	432.0	3	104	Callosobruchus	maculatus	half-sib	half-sib	insect	WB035.G	WB035.S	BS
./data_GP_compare/WB035/S1_n=2207_P.csv	Messina et al.	2013	WB035	S	P	2207	1319.5	2204.0	3	103	Callosobruchus	maculatus	half-sib	half-sib	insect	WB035.P	WB035.S	BS
./data_GP_compare/WB038/B_n=135_G.csv	Rauter&Moore	2002	WB038	B	G	135	387.0	135.0	4	105	Nicrophorus	pustulatus	half-sib	half-sib	insect	WB038.G	WB038.B	BS
./data_GP_compare/WB038/B_n=639_P.csv	Rauter&Moore	2002	WB038	B	P	639	387.0	639.0	4	106	Nicrophorus	pustulatus	half-sib	half-sib	insect	WB038.P	WB038.B	BS
./data_GP_compare/WB038/S1_n=135_G.csv	Rauter&Moore	2002	WB038	S	G	135	387.0	135.0	4	107	Nicrophorus	pustulatus	half-sib	half-sib	insect	WB038.G	WB038.S	BS
./data_GP_compare/WB038/S1_n=636_P.csv	Rauter&Moore	2002	WB038	S	P	636	387.0	639.0	4	108	Nicrophorus	pustulatus	half-sib	half-sib	insect	WB038.P	WB038.S	BS
./data_GP_compare/WB045/A_n=39_G.csv	Simons & Roff	1996	WB045	A	G	39	209.5	39.0	7	110	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045.G	WB045.A	AN
./data_GP_compare/WB045/A_n=380_P.csv	Simons & Roff	1996	WB045	A	P	380	209.5	371.0	7	109	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045.P	WB045.A	AN
./data_GP_compare/WB045/N1_n=39_G.csv	Simons & Roff	1996	WB045	N	G	39	200.5	39.0	7	112	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045.G	WB045.N	AN
./data_GP_compare/WB045/N1_n=362_P.csv	Simons & Roff	1996	WB045	N	P	362	200.5	371.0	7	111	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045.P	WB045.N	AN
./data_GP_compare/WB045a/A_n=39_G.csv	Simons & Roff	1996	WB045a	A	G	39	215.0	39.0	6	113	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045a.G	WB045a.A	AN
./data_GP_compare/WB045a/A_n=391_P.csv	Simons & Roff	1996	WB045a	A	P	391	215.0	381.0	6	114	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045a.P	WB045a.A	AN
./data_GP_compare/WB045a/N1_n=39_G.csv	Simons & Roff	1996	WB045a	N	G	39	205.0	39.0	6	116	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045a.G	WB045a.N	AN
./data_GP_compare/WB045a/N1_n=371_P.csv	Simons & Roff	1996	WB045a	N	P	371	205.0	381.0	6	115	Gryllus	pennsylvanicus	full-sib	full-sib	insect	WB045a.P	WB045a.N	AN
./data_GP_compare/WB053/A_n=24_G.csv	Via & Conner	1995	WB053	A	G	24	252.0	24.0	2	117	Tribolium	castaneum	half-sib	half-sib	insect	WB053.G	WB053.A	AN
./data_GP_compare/WB053/A_n=480_P.csv	Via & Conner	1995	WB053	A	P	480	252.0	480.0	2	118	Tribolium	castaneum	half-sib	half-sib	insect	WB053.P	WB053.A	AN
./data_GP_compare/WB053/N1_n=24_G.csv	Via & Conner	1995	WB053	N	G	24	252.0	24.0	2	119	Tribolium	castaneum	half-sib	half-sib	insect	WB053.G	WB053.N	AN
./data_GP_compare/WB053/N1_n=480_P.csv	Via & Conner	1995	WB053	N	P	480	252.0	480.0	2	120	Tribolium	castaneum	half-sib	half-sib	insect	WB053.P	WB053.N	AN
./data_GP_compare/WB053_3/A_n=24_G.csv	Via & Conner	1995	WB053_3	A	G	24	252.0	24.0	2	121	Tribolium	castaneum	half-sib	half-sib	insect	WB053_3.G	WB053_3.A	AN
./data_GP_compare/WB053_3/A_n=480_P.csv	Via & Conner	1995	WB053_3	A	P	480	252.0	480.0	2	122	Tribolium	castaneum	half-sib	half-sib	insect	WB053_3.P	WB053_3.A	AN
./data_GP_compare/WB053_3/N3_n=24_G.csv	Via & Conner	1995	WB053_3	N	G	24	252.0	24.0	2	123	Tribolium	castaneum	half-sib	half-sib	insect	WB053_3.G	WB053_3.N	AN
./data_GP_compare/WB053_3/N3_n=480_P.csv	Via & Conner	1995	WB053_3	N	P	480	252.0	480.0	2	124	Tribolium	castaneum	half-sib	half-sib	insect	WB053_3.P	WB053_3.N	AN
./data_GP_compare/WB053b/A_n=80_G.csv	Via & Conner	1995	WB053b	A	G	80	280.0	80.0	2	126	Tribolium	castaneum	half-sib	half-sib	insect	WB053b.G	WB053b.A	AN
./data_GP_compare/WB053b/A_n=480_P.csv	Via & Conner	1995	WB053b	A	P	480	280.0	480.0	2	125	Tribolium	castaneum	half-sib	half-sib	insect	WB053b.P	WB053b.A	AN
./data_GP_compare/WB053b/N1_n=80_G.csv	Via & Conner	1995	WB053b	N	G	80	280.0	80.0	2	128	Tribolium	castaneum	half-sib	half-sib	insect	WB053b.G	WB053b.N	AN
./data_GP_compare/WB053b/N1_n=480_P.csv	Via & Conner	1995	WB053b	N	P	480	280.0	480.0	2	127	Tribolium	castaneum	half-sib	half-sib	insect	WB053b.P	WB053b.N	AN
./data_GP_compare/WB053b_3/A_n=80_G.csv	Via & Conner	1995	WB053b_3	A	G	80	280.0	80.0	2	130	Tribolium	castaneum	half-sib	half-sib	insect	WB053b_3.G	WB053b_3.A	AN
./data_GP_compare/WB053b_3/A_n=480_P.csv	Via & Conner	1995	WB053b_3	A	P	480	280.0	480.0	2	129	Tribolium	castaneum	half-sib	half-sib	insect	WB053b_3.P	WB053b_3.A	AN
./data_GP_compare/WB053b_3/N3_n=80_G.csv	Via & Conner	1995	WB053b_3	N	G	80	280.0	80.0	2	132	Tribolium	castaneum	half-sib	half-sib	insect	WB053b_3.G	WB053b_3.N	AN
./data_GP_compare/WB053b_3/N3_n=480_P.csv	Via & Conner	1995	WB053b_3	N	P	480	280.0	480.0	2	131	Tribolium	castaneum	half-sib	half-sib	insect	WB053b_3.P	WB053b_3.N	AN
./data_GP_compare/WB057/A_n=22_G.csv	Grill et al.	1997	WB057	A	G	22	149.5	24.0	4	133	Harmonia	axyridis	full-sib	full-sib	insect	WB057.G	WB057.A	AN
./data_GP_compare/WB057/A_n=277_P.csv	Grill et al.	1997	WB057	A	P	277	149.5	277.0	4	134	Harmonia	axyridis	full-sib	full-sib	insect	WB057.P	WB057.A	AN
./data_GP_compare/WB057/N1_n=26_G.csv	Grill et al.	1997	WB057	N	G	26	151.5	24.0	4	135	Harmonia	axyridis	full-sib	full-sib	insect	WB057.G	WB057.N	AN
./data_GP_compare/WB057/N1_n=277_P.csv	Grill et al.	1997	WB057	N	P	277	151.5	277.0	4	136	Harmonia	axyridis	full-sib	full-sib	insect	WB057.P	WB057.N	AN

file – Described the path to the filename for the specific matrix that has been extracted from the paper.

authors – First name of author for the paper matrices were extracted from.

year – Year study was published

study – Study identifier

env – Environment traits were measured. B = Benign; S = stressful; A = non-novel; N = novel. Note that B and S are really just more detailed descriptors of the environments being stressful based on reading of the paper and the authors’ understanding of the system.

type – The type of matrix.G = Genetic variance/covariance matrix; P Phenotypic variance/covariance matrix

n – The average sample size. For P this is the total number of individuals, if G this is the number of families.

avg_n_env – Is the average sample size across environments

avg_n_type – Is the average sample size across the G and P matrices

trait_num – The number of traits in the matrices

col – Column indicator which matches names of matrices with meta-data.

genus Genus name

species Binomial species name

design_raw – The raw breeding design described in the paper. Can be half-sib, full-sib or RIL (Recombinant Inbred Lines). These are recoded so that RIL is just part of a half-sib design

design – Simplified breeding design

group – General taxonomic group studied

name.type – Interaction/ composite variable of study number and matrix type.

name.env – Interaction/ composite variable of study number and environment type.

env.comp – Environmental comparison moderator; BS = benign/stressful comparisons; AN = non-novel/novel comparison.

1.1.4 – Importing the Matrices

Import matrices and relevant data.

         mats <- import_mat(matrices)
        dscrp <- import_mat(descripDat)

Scale the data to prevent disproportionately large effects of the traits with large values. Not scaling will likely lead to similar results, but this is a safer option.

  standard <- list()
  for (i in 1:length(dscrp)){
    # standardize descriptive data
    # calculate the average mean for each trait
    dscrp[[i]]$stanfac <- ave(dscrp[[i]]$mean,dscrp[[i]]$trait)
    #standardize mean, sd, se, Va
    dscrp[[i]]$mean <- dscrp[[i]]$mean/dscrp[[i]]$stanfac
    dscrp[[i]]$sd <- dscrp[[i]]$sd/dscrp[[i]]$stanfac
    dscrp[[i]]$se <- dscrp[[i]]$se/dscrp[[i]]$stanfac
    dscrp[[i]]$Va <- dscrp[[i]]$Va/(dscrp[[i]]$stanfac^2)
    # save the standardization factor in to another list
    standard[[i]] <- data.frame(trait=unique(dscrp[[i]]$trait),stanfac=dscrp[[i]]$stanfac[which(!duplicated(dscrp[[i]]$trait))]) 
    # remove spaces and dashes from the trait names to be consistent with the names in the covariance-matrices.
    standard[[i]]$trait <- gsub(" ","\\.", as.character(standard[[i]]$trait))
    standard[[i]]$trait <- gsub("-","\\.", standard[[i]]$trait)
    rownames(standard[[i]]) <- standard[[i]]$trait
    }
  
  id_d <- sapply(names(dscrp),function(x){strsplit(x,"/")[[1]][3]})
  id_m <- sapply(names(mats),function(x){strsplit(x,"/")[[1]][3]})
  for (i in 1:length(mats)){
    # select the right standarization factors
    j <- which(id_d == id_m[[i]])
    std <- standard[[j]]
    # remove traits not represented in the covar-matrix
    std <- std[which(std$trait %in% names(mats[[i]])),]
    # order traits as in the covar-matrix
    std <- std[names(mats[[i]]),]
    # check whether the standardization factors line up
    if(sum(names(mats[[i]]) != std$trait)>0){stop("There is a inconsistency between trait names of a matrix and a description file")}
    mats[[i]] <- mats[[i]]/(std$stanfac %*% t(std$stanfac))
    }
  
  # remove studies without means which therefore cannot be standardized 
  rm.ids <- which(is.na(lapply(mats,sum)))
  mats <- mats[-rm.ids]
  dscrp <- dscrp[-which(id_d %in% id_m[rm.ids])]

Extract the sample size and get a sense of the number of traits from the filenames of each list.

                n <- match_n(names(mats))                       
        trait_num <- unlist(lapply(mats, function(x) nrow(x)))

1.2 – Simulations using G and P matrices with given sample sizes

Now that we have loaded all the data we can run the simulations described in the paper. This is done iteratively for each matrix and for each effect size. This generally does not take too long, even for 5000 simulations.

        sims = 5000

        system.time(
            # Proportion of variance of the dominant eigen value
        tota_var <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "totalVar", cor = FALSE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )

##    user  system elapsed 
## 407.424   9.352 420.286

        system.time(
            # Proportion of variance of the dominant eigen value
        var_along_logit <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "varAlong", cor = FALSE, logit = TRUE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )

##    user  system elapsed 
## 554.818   9.704 568.664

        system.time(
            # raw simulated matrices 
        matrices <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "raw_matrix", cor = FALSE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )

##    user  system elapsed 
## 200.958   7.285 208.637

1.3 – Generating Effect Sizes and Sampling Variances for Meta-analysis

Now that we have got all the simulations of the G and P matrices for the three different effect sizes of interest, we can extract from these distributions the average effect size and corresponding sampling error (i.e., the standard deviation). Note that, to generate some of the effect sizes we need to rely on the raw simulated matrices (i.e., the matrices object)

First, lets create the datasets for the first two effect sizes as they are already done for us. We will plot these data to to have a look and make sure nothing crazy is happening. We will also calculate the weights for the effect sizes and add these into the data frame.

First for standardised mean difference (SMD) in variance:

# Extract the relevant columns from the metadata file
     meta_sub_env <- meta[,c("authors", "year", "study",  "env", "avg_n_env", "trait_num", "genus", "species", "design", "group", "env.comp")]

    meta_sub_type <- meta[,c("authors", "year", "study", "type",  "avg_n_type", "trait_num", "genus", "species", "design", "group", "env.comp")]

We can have a look at the simulated distributions for each study graphically and print out a list used for comparisons to make sure things are working as expected.

    create_data(tota_var, 
                type = "SMD", 
                environment = "across", 
                plot = "yes")

Figure 1.3.1 – Effect size distributions of the standardised mean differences in total variance between matrices.

Now that this looks good we can create the data, add important information to the data frame, and save the datafile so that we no longer need to re-run these simulations each time we want to re-run analyses.

# Final data 
# Standardized mean difference in variance between G-G and P-P across environments
    
               SMD <- create_data(tota_var, 
                                  type = "SMD", 
                                  environment = "across", 
                                  plot = "no")
    
          SMD_data <- combine_meta_data(SMD, 
                                        meta_sub_type, 
                                        environment = "across")

    SMD_data$err   <- 1:nrow(SMD_data)
     SMD_data$wi   <- 1/SMD_data$S_var
    
    write.csv(SMD_data, file = "./data_GP_compare/SMD_data_standardised.csv")

Directionality of Effects: It’s important to remember here that positive values mean that the total variance in \(G_{novel}\) or \(P_{novel}\) is larger than \(G_{non-novel}\) or \(P_{non-novel}\), where as negative values mean that the total variance in \(G_{novel}\) or \(P_{novel}\) is smaller than \(G_{non-novel}\) or \(P_{non-novel}\).

Second, we can get the logit transformed proportion of variance along \(g_{max}\). This is to make the differences at the extremes (i.e., 0 and 1) less crazy compared to other values as it stretches out the ends of the distribution and normalises things as best as possible. Again, we can plot the simulations out to make sure that they look good before creating the data.

             create_data(var_along_logit, 
                        type = "raw_diff", 
                        environment = "across", 
                        plot = "yes")

Figure 1.3.2 – Effect size distributions of the logit transformed proportion of variance along \(g/p_{max}\) between matrices.

Once we are satisfied, we can create the data frame so we don’t need to run the simulations again.

          var_along_logit_d <- create_data(var_along_logit, 
                                           type = "raw_diff", 
                                           environment = "across", 
                                           plot = "no")
  
        var_along_logit_data <- combine_meta_data(var_along_logit_d, 
                                                   meta_sub_type, 
                                                   environment = "across")

    var_along_logit_data$err <- 1:nrow(var_along_logit_data)
     var_along_logit_data$wi <- 1/var_along_logit_data$S_var

    write.csv(var_along_logit_data, file = "./data_GP_compare/var_along_logit_data_standardised.csv")

Directionality of Effects: It’s important to remember here that positive values mean that the proportion of variance along \(gmax_{novel}\) or \(pmax_{novel}\) is larger than \(gmax_{non-novel}\) or \(pmax_{non-novel}\), where as negative values mean that the total variance in \(gmax_{novel}\) or \(pmax_{novel}\) is smaller than \(gmax_{non-novel}\) or \(pmax_{non-novel}\). We also need to remember that these are logit transformed.

Last, we can then calculate the angle between G-G across environments and G-P within each of the environments.

Plot to have a look at the effect size distribution across studies:

            create_data(matrices, 
                        type = "angle", 
                        environment = "across", 
                        plot = "yes")

Figure 1.3.3 – Effect size distributions of the logit transformed angle between \(g_{max}\) across environments.

The distributions are not as clean as the others, but not terrible and we can create the data for comparisons of the same matrix (e.g., \(G_{novel}\) compared to \(G_{non-novel}\)) across environments.

# Angle between P-P and G-G across environments for each study.
    
             angle_across <- create_data(matrices, 
                                            type = "angle", 
                                            environment = "across", 
                                            plot = "no")
    
        angle_across_data <- combine_meta_data(angle_across, 
                                                meta_sub_type, 
                                                environment = "across")

    angle_across_data$err <- 1:nrow(angle_across_data)
     angle_across_data$wi <- 1/angle_across_data$S_var

    write.csv(angle_across_data, file = "./data_GP_compare/angle_across_data_standardised.csv")

Directionality of Effects: Angles are logit transformed by log transformation of the angle itself and dividing this by 90 degrees. Positive effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(gmax_{novel}\) or \(pmax_{non-novel}\) and \(pmax_{novel}\) is approaching 90 degrees (i.e., values around 7 = 90 degrees). In contrast, negative effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(gmax_{novel}\) or \(pmax_{non-novel}\) and \(pmax_{novel}\) is approaching 0 degrees (i.e., values around -7 = 0 degrees).

We can then do the same for comparisons of different matrices within each environment to understand how G and P align, for example (i.e., \(G_{novel}\) compared to \(P_{novel}\))

        create_data(matrices, type = "angle", environment = "within", plot = "yes")

Figure 1.3.4 – Effect size distributions of the logit transformed angle between \(g_{max}\) and \(p_{max}\) within environments.

Once we are happy with this we can create the data frame so we again don’t need to re-run these simulations again.

    # Angle between P and G within an environment for each study.

             angle_within <- create_data(matrices, type = "angle", environment = "within", plot = "no")
    
        angle_within_data <- combine_meta_data(angle_within, meta_sub_env, environment = "within")

    angle_within_data$err <- 1:nrow(angle_within_data)
     angle_within_data$wi <- 1/angle_within_data$S_var

    write.csv(angle_within_data, file = "./data_GP_compare/angle_within_data_standardised.csv")

Directionality of Effects: Angles are logit transformed by log transformation of the angle itself and dividing this by 90 degrees. Positive effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(pmax_{non-novel}\) or \(gmax_{novel}\) and \(pmax_{novel}\) is approaching 90 degrees (i.e., values around 7 = 90 degrees). In contrast, negative effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(pmax_{non-novel}\) or \(gmax_{novel}\) and \(pmax_{novel}\) is approaching 0 degrees (i.e., values around -7 = 0 degrees). Values in between are varying angles.

2.0 – Simulation 2: Alignment of Plastic Responses with G_novel and G_non-novel

Now that we have simulated all the data that is relevant for comparisons between G and P we can begin to ask questions about the evolvability of plastic responses. Since we only need the mean trait changes for traits across environments to generate a phenotype response vector we can actually increase the total sample size of studies included because we don’t need to rely on getting the entire P matrix for simulations.

2.1 – Importing the data

To import the data we have separated out all of the matrices for each of the two simulations and put these in separate folders. This just makes it a little more clear what data was used for each of the simulations described.

      matrices_evolve <- files(dir = "./data_evolvability/", type = "matrix", file_pattern = "[WRs][BR]")
    descripDat_evolve <- files(dir = "./data_evolvability/", type = "D", file_pattern = "[WRs][BR]")

We also need to bring in some new meta-data that includes the additional studies that are added. One can view the data below.

    meta_evolve <- read.csv("./data_evolvability/metadata_final.csv", stringsAsFactor = FALSE)
    kable(meta_evolve) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "200px")

file	authors	year	study	env	type	n	trait_num	col	genus	species	design_raw	design	group	env.comp	name.type	name.env
./data_evolvability/RR062/A_n=10_G.csv	Beacham	1988	RR062	A	G	10	3	1	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062.G	RR062.A
./data_evolvability/RR062/A_n=100_P.csv	Beacham	1988	RR062	A	P	100	3	2	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062.P	RR062.A
./data_evolvability/RR062/N1_n=10_G.csv	Beacham	1988	RR062	N	G	10	3	3	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062.G	RR062.N
./data_evolvability/RR062/N1_n=74_P.csv	Beacham	1988	RR062	N	P	74	3	4	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062.P	RR062.N
./data_evolvability/RR062_b/A_n=10_G.csv	Beacham	1988	RR062_b	A	G	10	3	5	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062_b.G	RR062_b.A
./data_evolvability/RR062_b/A_n=100_P.csv	Beacham	1988	RR062_b	A	P	100	3	6	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062_b.P	RR062_b.A
./data_evolvability/RR062_b/N2_n=10_G.csv	Beacham	1988	RR062_b	N	G	10	3	7	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062_b.G	RR062_b.N
./data_evolvability/RR062_b/N2_n=100_P.csv	Beacham	1988	RR062_b	N	P	100	3	8	Oncorhynchus	gorbuscha	half-sib	half-sib	fish	AN	RR062_b.P	RR062_b.N
./data_evolvability/RR062b/A_n=10_G.csv	Beacham	1988	RR062b	A	G	10	3	9	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b.G	RR062b.A
./data_evolvability/RR062b/A_n=100_P.csv	Beacham	1988	RR062b	A	P	100	3	10	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b.P	RR062b.A
./data_evolvability/RR062b/N1_n=10_G.csv	Beacham	1988	RR062b	N	G	10	3	11	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b.G	RR062b.N
./data_evolvability/RR062b/N1_n=98_P.csv	Beacham	1988	RR062b	N	P	98	3	12	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b.P	RR062b.N
./data_evolvability/RR062b_b/A_n=10_G.csv	Beacham	1988	RR062b_b	A	G	10	3	13	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b_b.G	RR062b_b.A
./data_evolvability/RR062b_b/A_n=100_P.csv	Beacham	1988	RR062b_b	A	P	100	3	14	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b_b.P	RR062b_b.A
./data_evolvability/RR062b_b/N2_n=10_G.csv	Beacham	1988	RR062b_b	N	G	10	3	15	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b_b.G	RR062b_b.N
./data_evolvability/RR062b_b/N2_n=98_P.csv	Beacham	1988	RR062b_b	N	P	98	3	16	Oncorhynchus	keta	half-sib	half-sib	fish	AN	RR062b_b.P	RR062b_b.N
./data_evolvability/RR073/B_n=68_G.csv	Evans et al.	2015	RR073	B	G	68	8	18	Poecilia	reticulata	half-sib	half-sib	fish	BS	RR073.G	RR073.B
./data_evolvability/RR073/B_n=161_P.csv	Evans et al.	2015	RR073	B	P	161	8	17	Poecilia	reticulata	half-sib	half-sib	fish	BS	RR073.P	RR073.B
./data_evolvability/RR073/S1_n=64_G.csv	Evans et al.	2015	RR073	S	G	64	8	20	Poecilia	reticulata	half-sib	half-sib	fish	BS	RR073.G	RR073.S
./data_evolvability/RR073/S1_n=146_P.csv	Evans et al.	2015	RR073	S	P	146	8	19	Poecilia	reticulata	half-sib	half-sib	fish	BS	RR073.P	RR073.S
./data_evolvability/RR076/B_n=69_G.csv	Guan et al.	2016	RR076	B	G	69	3	22	Scophthalmus	maximus	half-sib	half-sib	fish	BS	RR076.G	RR076.B
./data_evolvability/RR076/B_n=2125_P.csv	Guan et al.	2016	RR076	B	P	2125	3	21	Scophthalmus	maximus	half-sib	half-sib	fish	BS	RR076.P	RR076.B
./data_evolvability/RR076/S1_n=69_G.csv	Guan et al.	2016	RR076	S	G	69	3	24	Scophthalmus	maximus	half-sib	half-sib	fish	BS	RR076.G	RR076.S
./data_evolvability/RR076/S1_n=2925_P.csv	Guan et al.	2016	RR076	S	P	2925	3	23	Scophthalmus	maximus	half-sib	half-sib	fish	BS	RR076.P	RR076.S
./data_evolvability/RR087a/A_n=100_G.csv	Sae-Lim et al.	2013	RR087a	A	G	100	3	25	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087a.G	RR087a.A
./data_evolvability/RR087a/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087a	A	P	2364	3	26	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087a.P	RR087a.A
./data_evolvability/RR087a/N1_n=100_G.csv	Sae-Lim et al.	2013	RR087a	N	G	100	3	27	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087a.G	RR087a.N
./data_evolvability/RR087a/N1_n=1891_P.csv	Sae-Lim et al.	2013	RR087a	N	P	1891	3	28	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087a.P	RR087a.N
./data_evolvability/RR087b/A_n=100_G.csv	Sae-Lim et al.	2013	RR087b	A	G	100	3	29	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087b.G	RR087b.A
./data_evolvability/RR087b/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087b	A	P	2364	3	30	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087b.P	RR087b.A
./data_evolvability/RR087b/N2_n=100_G.csv	Sae-Lim et al.	2013	RR087b	N	G	100	3	31	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087b.G	RR087b.N
./data_evolvability/RR087b/N2_n=2790_P.csv	Sae-Lim et al.	2013	RR087b	N	P	2790	3	32	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087b.P	RR087b.N
./data_evolvability/RR087c/A_n=100_G.csv	Sae-Lim et al.	2013	RR087c	A	G	100	3	33	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087c.G	RR087c.A
./data_evolvability/RR087c/A_n=2364_P.csv	Sae-Lim et al.	2013	RR087c	A	P	2364	3	34	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087c.P	RR087c.A
./data_evolvability/RR087c/N3_n=100_G.csv	Sae-Lim et al.	2013	RR087c	N	G	100	3	35	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087c.G	RR087c.N
./data_evolvability/RR087c/N3_n=1818_P.csv	Sae-Lim et al.	2013	RR087c	N	P	1818	3	36	Oncorhynchus	mykiss	full-sib	full-sib	fish	AN	RR087c.P	RR087c.N
./data_evolvability/sR097/B_n=26_G.csv	Kasule	1991	sR097	B	G	26	2	37	Dysdercus	fasciatus	half-sib	half-sib	insect	BS	sR097.G	sR097.B
./data_evolvability/sR097/B_n=312_P.csv	Kasule	1991	sR097	B	P	312	2	38	Dysdercus	fasciatus	half-sib	half-sib	insect	BS	sR097.P	sR097.B
./data_evolvability/sR097/S1_n=26_G.csv	Kasule	1991	sR097	S	G	26	2	39	Dysdercus	fasciatus	half-sib	half-sib	insect	BS	sR097.G	sR097.S
./data_evolvability/sR097/S1_n=312_P.csv	Kasule	1991	sR097	S	P	312	2	40	Dysdercus	fasciatus	half-sib	half-sib	insect	BS	sR097.P	sR097.S
./data_evolvability/sR110b/B_n=37_G.csv	Blanckenhorn & Heyland	2004	sR110b	B	G	37	2	42	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110b.G	sR110b.B
./data_evolvability/sR110b/B_n=128_P.csv	Blanckenhorn & Heyland	2004	sR110b	B	P	128	2	41	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110b.P	sR110b.B
./data_evolvability/sR110b/S1_n=27_G.csv	Blanckenhorn & Heyland	2004	sR110b	S	G	27	2	43	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110b.G	sR110b.S
./data_evolvability/sR110b/S1_n=98_P.csv	Blanckenhorn & Heyland	2004	sR110b	S	P	98	2	44	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110b.P	sR110b.S
./data_evolvability/sR110c/B_n=40_G.csv	Blanckenhorn & Heyland	2004	sR110c	B	G	40	2	46	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110c.G	sR110c.B
./data_evolvability/sR110c/B_n=171_P.csv	Blanckenhorn & Heyland	2004	sR110c	B	P	171	2	45	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110c.P	sR110c.B
./data_evolvability/sR110c/S1_n=33_G.csv	Blanckenhorn & Heyland	2004	sR110c	S	G	33	2	48	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110c.G	sR110c.S
./data_evolvability/sR110c/S1_n=129_P.csv	Blanckenhorn & Heyland	2004	sR110c	S	P	129	2	47	Scathophaga	stercoraria	half-sib	half-sib	insect	BS	sR110c.P	sR110c.S
./data_evolvability/sR115/A_n=350_G.csv	Conner et al.	2003	sR115	A	G	350	6	50	Raphanus	raphanistrum	half-sib	half-sib	plant	AN	sR115.G	sR115.A
./data_evolvability/sR115/A_n=1133_P.csv	Conner et al.	2003	sR115	A	P	1133	6	49	Raphanus	raphanistrum	half-sib	half-sib	plant	AN	sR115.P	sR115.A
./data_evolvability/sR115/N_n=593_G.csv	Conner et al.	2003	sR115	N	G	593	6	51	Raphanus	raphanistrum	half-sib	half-sib	plant	AN	sR115.G	sR115.N
./data_evolvability/sR115/N_n=888_P.csv	Conner et al.	2003	sR115	N	P	888	6	52	Raphanus	raphanistrum	half-sib	half-sib	plant	AN	sR115.P	sR115.N
./data_evolvability/sR116/A_n=366_G.csv	Delcourt et al.	2010	sR116	A	G	366	8	54	Drosophila	serrata	half-sib	half-sib	insect	AN	sR116.G	sR116.A
./data_evolvability/sR116/A_n=3315_P.csv	Delcourt et al.	2010	sR116	A	P	3315	8	53	Drosophila	serrata	half-sib	half-sib	insect	AN	sR116.P	sR116.A
./data_evolvability/sR116/N1_n=365_G.csv	Delcourt et al.	2010	sR116	N	G	365	8	56	Drosophila	serrata	half-sib	half-sib	insect	AN	sR116.G	sR116.N
./data_evolvability/sR116/N1_n=3283_P.csv	Delcourt et al.	2010	sR116	N	P	3283	8	55	Drosophila	serrata	half-sib	half-sib	insect	AN	sR116.P	sR116.N
./data_evolvability/sR145/B_n=52_G.csv	Tucic et al.	1991	sR145	B	G	52	8	58	Acanthoscelides	obtectus	full-sib	full-sib	insect	BS	sR145.G	sR145.B
./data_evolvability/sR145/B_n=150_P.csv	Tucic et al.	1991	sR145	B	P	150	8	57	Acanthoscelides	obtectus	full-sib	full-sib	insect	BS	sR145.P	sR145.B
./data_evolvability/sR145/S1_n=43_G.csv	Tucic et al.	1991	sR145	S	G	43	8	60	Acanthoscelides	obtectus	full-sib	full-sib	insect	BS	sR145.G	sR145.S
./data_evolvability/sR145/S1_n=124_P.csv	Tucic et al.	1991	sR145	S	P	124	8	59	Acanthoscelides	obtectus	full-sib	full-sib	insect	BS	sR145.P	sR145.S
./data_evolvability/WB001/B_n=30_G.csv	Auld	2010	WB001	B	G	30	3	62	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001.G	WB001.B
./data_evolvability/WB001/B_n=210_P.csv	Auld	2010	WB001	B	P	210	3	61	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001.P	WB001.B
./data_evolvability/WB001/S1_n=30_G.csv	Auld	2010	WB001	S	G	30	3	64	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001.G	WB001.S
./data_evolvability/WB001/S1_n=150_P.csv	Auld	2010	WB001	S	P	150	3	63	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001.P	WB001.S
./data_evolvability/WB001_b/B_n=30_G.csv	Auld	2010	WB001_b	B	G	30	3	66	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_b.G	WB001_b.B
./data_evolvability/WB001_b/B_n=210_P.csv	Auld	2010	WB001_b	B	P	210	3	65	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_b.P	WB001_b.B
./data_evolvability/WB001_b/S2_n=30_G.csv	Auld	2010	WB001_b	S	G	30	3	68	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_b.G	WB001_b.S
./data_evolvability/WB001_b/S2_n=210_P.csv	Auld	2010	WB001_b	S	P	210	3	67	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_b.P	WB001_b.S
./data_evolvability/WB001_c/B_n=30_G.csv	Auld	2010	WB001_c	B	G	30	3	70	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_c.G	WB001_c.B
./data_evolvability/WB001_c/B_n=210_P.csv	Auld	2010	WB001_c	B	P	210	3	69	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_c.P	WB001_c.B
./data_evolvability/WB001_c/S3_n=30_G.csv	Auld	2010	WB001_c	S	G	30	3	72	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_c.G	WB001_c.S
./data_evolvability/WB001_c/S3_n=150_P.csv	Auld	2010	WB001_c	S	P	150	3	71	Physa	acuta	full-sib	full-sib	mollusc	BS	WB001_c.P	WB001_c.S
./data_evolvability/WB002a/A_n=39_G.csv	Begin & Roff	2001	WB002a	A	G	39	5	73	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB002a.G	WB002a.A
./data_evolvability/WB002a/A_n=705_P.csv	Begin & Roff	2001	WB002a	A	P	705	5	74	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB002a.P	WB002a.A
./data_evolvability/WB002a/N1_n=39_G.csv	Begin & Roff	2001	WB002a	N	G	39	5	75	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB002a.G	WB002a.N
./data_evolvability/WB002a/N1_n=505_P.csv	Begin & Roff	2001	WB002a	N	P	505	5	76	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB002a.P	WB002a.N
./data_evolvability/WB003/A_n=56_G.csv	Begin et al.	2004	WB003	A	G	56	5	78	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003.G	WB003.A
./data_evolvability/WB003/A_n=533_P.csv	Begin et al.	2004	WB003	A	P	533	5	77	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003.P	WB003.A
./data_evolvability/WB003/N1_n=62_G.csv	Begin et al.	2004	WB003	N	G	62	5	79	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003.G	WB003.N
./data_evolvability/WB003/N1_n=862_P.csv	Begin et al.	2004	WB003	N	P	862	5	80	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003.P	WB003.N
./data_evolvability/WB003_2/A_n=56_G.csv	Begin et al.	2004	WB003_2	A	G	56	5	82	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003_2.G	WB003_2.A
./data_evolvability/WB003_2/A_n=533_P.csv	Begin et al.	2004	WB003_2	A	P	533	5	81	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003_2.P	WB003_2.A
./data_evolvability/WB003_2/N2_n=60_G.csv	Begin et al.	2004	WB003_2	N	G	60	5	83	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003_2.G	WB003_2.N
./data_evolvability/WB003_2/N2_n=715_P.csv	Begin et al.	2004	WB003_2	N	P	715	5	84	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003_2.P	WB003_2.N
./data_evolvability/WB003a/A_n=45_G.csv	Begin et al.	2004	WB003a	A	G	45	5	86	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a.G	WB003a.A
./data_evolvability/WB003a/A_n=262_P.csv	Begin et al.	2004	WB003a	A	P	262	5	85	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a.P	WB003a.A
./data_evolvability/WB003a/N1_n=41_G.csv	Begin et al.	2004	WB003a	N	G	41	5	88	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a.G	WB003a.N
./data_evolvability/WB003a/N1_n=315_P.csv	Begin et al.	2004	WB003a	N	P	315	5	87	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a.P	WB003a.N
./data_evolvability/WB003a_2/A_n=45_G.csv	Begin et al.	2004	WB003a_2	A	G	45	5	90	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a_2.G	WB003a_2.A
./data_evolvability/WB003a_2/A_n=262_P.csv	Begin et al.	2004	WB003a_2	A	P	262	5	89	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a_2.P	WB003a_2.A
./data_evolvability/WB003a_2/N2_n=32_G.csv	Begin et al.	2004	WB003a_2	N	G	32	5	92	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a_2.G	WB003a_2.N
./data_evolvability/WB003a_2/N2_n=133_P.csv	Begin et al.	2004	WB003a_2	N	P	133	5	91	Gryllus	firmus	full-sib	full-sib	insect	AN	WB003a_2.P	WB003a_2.N
./data_evolvability/WB007/B_n=132_G.csv	Brock et al.	2010	WB007	B	G	132	11	93	Brassica	rapa	RIL	half-sib	plant	BS	WB007.G	WB007.B
./data_evolvability/WB007/B_n=601_P.csv	Brock et al.	2010	WB007	B	P	601	11	94	Brassica	rapa	RIL	half-sib	plant	BS	WB007.P	WB007.B
./data_evolvability/WB007/S1_n=131_G.csv	Brock et al.	2010	WB007	S	G	131	11	95	Brassica	rapa	RIL	half-sib	plant	BS	WB007.G	WB007.S
./data_evolvability/WB007/S1_n=524_P.csv	Brock et al.	2010	WB007	S	P	524	11	96	Brassica	rapa	RIL	half-sib	plant	BS	WB007.P	WB007.S
./data_evolvability/WB010/A_n=141_G.csv	Collins et al.	1999	WB010	A	G	141	4	97	Achroia	grisella	half-sib	half-sib	insect	AN	WB010.G	WB010.A
./data_evolvability/WB010/A_n=1540_P.csv	Collins et al.	1999	WB010	A	P	1540	4	98	Achroia	grisella	half-sib	half-sib	insect	AN	WB010.P	WB010.A
./data_evolvability/WB010/N1_n=62_G.csv	Collins et al.	1999	WB010	N	G	62	4	100	Achroia	grisella	half-sib	half-sib	insect	AN	WB010.G	WB010.N
./data_evolvability/WB010/N1_n=443_P.csv	Collins et al.	1999	WB010	N	P	443	4	99	Achroia	grisella	half-sib	half-sib	insect	AN	WB010.P	WB010.N
./data_evolvability/WB011/B_n=404_G.csv	Czesak & Fox	2003	WB011	B	G	404	2	102	Stator	limbatus	half-sib	half-sib	insect	BS	WB011.G	WB011.B
./data_evolvability/WB011/B_n=2667_P.csv	Czesak & Fox	2003	WB011	B	P	2667	2	101	Stator	limbatus	half-sib	half-sib	insect	BS	WB011.P	WB011.B
./data_evolvability/WB011/S1_n=404_G.csv	Czesak & Fox	2003	WB011	S	G	404	2	104	Stator	limbatus	half-sib	half-sib	insect	BS	WB011.G	WB011.S
./data_evolvability/WB011/S1_n=2667_P.csv	Czesak & Fox	2003	WB011	S	P	2667	2	103	Stator	limbatus	half-sib	half-sib	insect	BS	WB011.P	WB011.S
./data_evolvability/WB014/A_n=1007_P.csv	Delcourt et al.	2009	WB014	A	P	1007	2	105	Drosophila	serrata	half-sib	half-sib	insect	AN	WB014.P	WB014.A
./data_evolvability/WB014/A_n=319_G.csv	Delcourt et al.	2009	WB014	A	G	319	2	106	Drosophila	serrata	half-sib	half-sib	insect	AN	WB014.G	WB014.A
./data_evolvability/WB014/N1_n=317_G.csv	Delcourt et al.	2009	WB014	N	G	317	2	107	Drosophila	serrata	half-sib	half-sib	insect	AN	WB014.G	WB014.N
./data_evolvability/WB014/N1_n=998_P.csv	Delcourt et al.	2009	WB014	N	P	998	2	108	Drosophila	serrata	half-sib	half-sib	insect	AN	WB014.P	WB014.N
./data_evolvability/WB016/B_n=27_G.csv	Engqvist	2007	WB016	B	G	27	2	110	Panorpa	cognata	full-sib	full-sib	insect	BS	WB016.G	WB016.B
./data_evolvability/WB016/B_n=140_P.csv	Engqvist	2007	WB016	B	P	140	2	109	Panorpa	cognata	full-sib	full-sib	insect	BS	WB016.P	WB016.B
./data_evolvability/WB016/S1_n=27_G.csv	Engqvist	2007	WB016	S	G	27	2	112	Panorpa	cognata	full-sib	full-sib	insect	BS	WB016.G	WB016.S
./data_evolvability/WB016/S1_n=140_P.csv	Engqvist	2007	WB016	S	P	140	2	111	Panorpa	cognata	full-sib	full-sib	insect	BS	WB016.P	WB016.S
./data_evolvability/WB022/A_n=40_G.csv	Guntrip et al.	1997	WB022	A	G	40	2	113	Callosobruchus	maculatus	half-sib	half-sib	insect	AN	WB022.G	WB022.A
./data_evolvability/WB022/A_n=531_P.csv	Guntrip et al.	1997	WB022	A	P	531	2	114	Callosobruchus	maculatus	half-sib	half-sib	insect	AN	WB022.P	WB022.A
./data_evolvability/WB022/N1_n=40_G.csv	Guntrip et al.	1997	WB022	N	G	40	2	115	Callosobruchus	maculatus	half-sib	half-sib	insect	AN	WB022.G	WB022.N
./data_evolvability/WB022/N1_n=505_P.csv	Guntrip et al.	1997	WB022	N	P	505	2	116	Callosobruchus	maculatus	half-sib	half-sib	insect	AN	WB022.P	WB022.N
./data_evolvability/WB025/A_n=65_G.csv	Holloway et al.	1990	WB025	A	G	65	4	117	Sitophilus	oryzae	half-sib	half-sib	insect	AN	WB025.G	WB025.A
./data_evolvability/WB025/A_n=823_P.csv	Holloway et al.	1990	WB025	A	P	823	4	118	Sitophilus	oryzae	half-sib	half-sib	insect	AN	WB025.P	WB025.A
./data_evolvability/WB025/N1_n=69_G.csv	Holloway et al.	1990	WB025	N	G	69	4	119	Sitophilus	oryzae	half-sib	half-sib	insect	AN	WB025.G	WB025.N
./data_evolvability/WB025/N1_n=778_P.csv	Holloway et al.	1990	WB025	N	P	778	4	120	Sitophilus	oryzae	half-sib	half-sib	insect	AN	WB025.P	WB025.N
./data_evolvability/WB028/B_n=60_G.csv	Klause & Morin	2001	WB028	B	G	60	2	121	Priophorus	pallipes	full-sib	full-sib	insect	BS	WB028.G	WB028.B
./data_evolvability/WB028/B_n=600_P.csv	Klause & Morin	2001	WB028	B	P	600	2	122	Priophorus	pallipes	full-sib	full-sib	insect	BS	WB028.P	WB028.B
./data_evolvability/WB028/S1_n=60_G.csv	Klause & Morin	2001	WB028	S	G	60	2	123	Priophorus	pallipes	full-sib	full-sib	insect	BS	WB028.G	WB028.S
./data_evolvability/WB028/S1_n=600_P.csv	Klause & Morin	2001	WB028	S	P	600	2	124	Priophorus	pallipes	full-sib	full-sib	insect	BS	WB028.P	WB028.S
./data_evolvability/WB029/B_n=63_G.csv	King et al.	2011	WB029	B	G	63	3	126	Gryllus	firmus	half-sib	half-sib	insect	BS	WB029.G	WB029.B
./data_evolvability/WB029/B_n=1662_P.csv	King et al.	2011	WB029	B	P	1662	3	125	Gryllus	firmus	half-sib	half-sib	insect	BS	WB029.P	WB029.B
./data_evolvability/WB029/S1_n=63_G.csv	King et al.	2011	WB029	S	G	63	3	128	Gryllus	firmus	half-sib	half-sib	insect	BS	WB029.G	WB029.S
./data_evolvability/WB029/S1_n=1728_P.csv	King et al.	2011	WB029	S	P	1728	3	127	Gryllus	firmus	half-sib	half-sib	insect	BS	WB029.P	WB029.S
./data_evolvability/WB033/A_n=60_G.csv	Lau et al.	2014	WB033	A	G	60	4	130	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033.G	WB033.A
./data_evolvability/WB033/A_n=526_P.csv	Lau et al.	2014	WB033	A	P	526	4	129	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033.P	WB033.A
./data_evolvability/WB033/N1_n=60_G.csv	Lau et al.	2014	WB033	N	G	60	4	132	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033.G	WB033.N
./data_evolvability/WB033/N1_n=500_P.csv	Lau et al.	2014	WB033	N	P	500	4	131	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033.P	WB033.N
./data_evolvability/WB033a/A_n=60_G.csv	Lau et al.	2014	WB033a	A	G	60	4	134	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033a.G	WB033a.A
./data_evolvability/WB033a/A_n=459_P.csv	Lau et al.	2014	WB033a	A	P	459	4	133	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033a.P	WB033a.A
./data_evolvability/WB033a/N1_n=60_G.csv	Lau et al.	2014	WB033a	N	G	60	4	136	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033a.G	WB033a.N
./data_evolvability/WB033a/N1_n=470_P.csv	Lau et al.	2014	WB033a	N	P	470	4	135	Arabidopsis	thaliana	RIL	half-sib	plant	AN	WB033a.P	WB033a.N
./data_evolvability/WB035/B_n=432_G.csv	Messina et al.	2013	WB035	B	G	432	3	138	Callosobruchus	maculatus	half-sib	half-sib	insect	BS	WB035.G	WB035.B
./data_evolvability/WB035/B_n=2201_P.csv	Messina et al.	2013	WB035	B	P	2201	3	137	Callosobruchus	maculatus	half-sib	half-sib	insect	BS	WB035.P	WB035.B
./data_evolvability/WB035/S1_n=432_G.csv	Messina et al.	2013	WB035	S	G	432	3	140	Callosobruchus	maculatus	half-sib	half-sib	insect	BS	WB035.G	WB035.S
./data_evolvability/WB035/S1_n=2207_P.csv	Messina et al.	2013	WB035	S	P	2207	3	139	Callosobruchus	maculatus	half-sib	half-sib	insect	BS	WB035.P	WB035.S
./data_evolvability/WB036/A_n=22_G.csv	Paccard et al.	2013	WB036	A	G	22	7	141	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036.G	WB036.A
./data_evolvability/WB036/A_n=58_P.csv	Paccard et al.	2013	WB036	A	P	58	7	142	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036.P	WB036.A
./data_evolvability/WB036/N1_n=22_G.csv	Paccard et al.	2013	WB036	N	G	22	7	143	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036.G	WB036.N
./data_evolvability/WB036/N1_n=63_P.csv	Paccard et al.	2013	WB036	N	P	63	7	144	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036.P	WB036.N
./data_evolvability/WB036A/A_n=22_G.csv	Paccard et al.	2013	WB036A	A	G	22	7	145	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036A.G	WB036A.A
./data_evolvability/WB036A/A_n=60_P.csv	Paccard et al.	2013	WB036A	A	P	60	7	146	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036A.P	WB036A.A
./data_evolvability/WB036A/N1_n=22_G.csv	Paccard et al.	2013	WB036A	N	G	22	7	147	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036A.G	WB036A.N
./data_evolvability/WB036A/N1_n=63_P.csv	Paccard et al.	2013	WB036A	N	P	63	7	148	Arabidopsis	lyrata	full-sib	full-sib	plant	AN	WB036A.P	WB036A.N
./data_evolvability/WB038/B_n=135_G.csv	Rauter&Moore	2002	WB038	B	G	135	4	149	Nicrophorus	pustulatus	half-sib	half-sib	insect	BS	WB038.G	WB038.B
./data_evolvability/WB038/B_n=639_P.csv	Rauter&Moore	2002	WB038	B	P	639	4	150	Nicrophorus	pustulatus	half-sib	half-sib	insect	BS	WB038.P	WB038.B
./data_evolvability/WB038/S1_n=135_G.csv	Rauter&Moore	2002	WB038	S	G	135	4	151	Nicrophorus	pustulatus	half-sib	half-sib	insect	BS	WB038.G	WB038.S
./data_evolvability/WB038/S1_n=636_P.csv	Rauter&Moore	2002	WB038	S	P	636	4	152	Nicrophorus	pustulatus	half-sib	half-sib	insect	BS	WB038.P	WB038.S
./data_evolvability/WB039/B_n=21_G.csv	Relyea	2005	WB039	B	G	21	9	153	Rana	sylvatica	half-sib	half-sib	amphibian	BS	WB039.G	WB039.B
./data_evolvability/WB039/B_n=2700_P.csv	Relyea	2005	WB039	B	P	2700	9	154	Rana	sylvatica	half-sib	half-sib	amphibian	BS	WB039.P	WB039.B
./data_evolvability/WB039/S1_n=21_G.csv	Relyea	2005	WB039	S	G	21	9	155	Rana	sylvatica	half-sib	half-sib	amphibian	BS	WB039.G	WB039.S
./data_evolvability/WB039/S1_n=2700_P.csv	Relyea	2005	WB039	S	P	2700	9	156	Rana	sylvatica	half-sib	half-sib	amphibian	BS	WB039.P	WB039.S
./data_evolvability/WB044/B_n=26_G.csv	Sherrard et al.	2009	WB044	B	G	26	10	158	Avena	barbata	RIL	half-sib	plant	BS	WB044.G	WB044.B
./data_evolvability/WB044/B_n=112_P.csv	Sherrard et al.	2009	WB044	B	P	112	10	157	Avena	barbata	RIL	half-sib	plant	BS	WB044.P	WB044.B
./data_evolvability/WB044/S1_n=26_G.csv	Sherrard et al.	2009	WB044	S	G	26	10	160	Avena	barbata	RIL	half-sib	plant	BS	WB044.G	WB044.S
./data_evolvability/WB044/S1_n=112_P.csv	Sherrard et al.	2009	WB044	S	P	112	10	159	Avena	barbata	RIL	half-sib	plant	BS	WB044.P	WB044.S
./data_evolvability/WB045/A_n=39_G.csv	Simons & Roff	1996	WB045	A	G	39	7	162	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045.G	WB045.A
./data_evolvability/WB045/A_n=380_P.csv	Simons & Roff	1996	WB045	A	P	380	7	161	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045.P	WB045.A
./data_evolvability/WB045/N1_n=39_G.csv	Simons & Roff	1996	WB045	N	G	39	7	164	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045.G	WB045.N
./data_evolvability/WB045/N1_n=362_P.csv	Simons & Roff	1996	WB045	N	P	362	7	163	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045.P	WB045.N
./data_evolvability/WB045a/A_n=39_G.csv	Simons & Roff	1996	WB045a	A	G	39	6	165	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045a.G	WB045a.A
./data_evolvability/WB045a/A_n=391_P.csv	Simons & Roff	1996	WB045a	A	P	391	6	166	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045a.P	WB045a.A
./data_evolvability/WB045a/N1_n=39_G.csv	Simons & Roff	1996	WB045a	N	G	39	6	168	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045a.G	WB045a.N
./data_evolvability/WB045a/N1_n=371_P.csv	Simons & Roff	1996	WB045a	N	P	371	6	167	Gryllus	pennsylvanicus	full-sib	full-sib	insect	AN	WB045a.P	WB045a.N
./data_evolvability/WB046/B_n=50_G.csv	Simonsen & Stinchcombe	2001	WB046	B	G	50	3	169	Ipomoea	hederacea	RIL	half-sib	plant	BS	WB046.G	WB046.B
./data_evolvability/WB046/B_n=600_P.csv	Simonsen & Stinchcombe	2001	WB046	B	P	600	3	170	Ipomoea	hederacea	RIL	half-sib	plant	BS	WB046.P	WB046.B
./data_evolvability/WB046/S1_n=50_G.csv	Simonsen & Stinchcombe	2001	WB046	S	G	50	3	171	Ipomoea	hederacea	RIL	half-sib	plant	BS	WB046.G	WB046.S
./data_evolvability/WB046/S1_n=600_P.csv	Simonsen & Stinchcombe	2001	WB046	S	P	600	3	172	Ipomoea	hederacea	RIL	half-sib	plant	BS	WB046.P	WB046.S
./data_evolvability/WB053/A_n=24_G.csv	Via & Conner	1995	WB053	A	G	24	2	173	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053.G	WB053.A
./data_evolvability/WB053/A_n=480_P.csv	Via & Conner	1995	WB053	A	P	480	2	174	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053.P	WB053.A
./data_evolvability/WB053/N1_n=24_G.csv	Via & Conner	1995	WB053	N	G	24	2	175	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053.G	WB053.N
./data_evolvability/WB053/N1_n=480_P.csv	Via & Conner	1995	WB053	N	P	480	2	176	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053.P	WB053.N
./data_evolvability/WB053_3/A_n=24_G.csv	Via & Conner	1995	WB053_3	A	G	24	2	177	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053_3.G	WB053_3.A
./data_evolvability/WB053_3/A_n=480_P.csv	Via & Conner	1995	WB053_3	A	P	480	2	178	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053_3.P	WB053_3.A
./data_evolvability/WB053_3/N3_n=24_G.csv	Via & Conner	1995	WB053_3	N	G	24	2	179	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053_3.G	WB053_3.N
./data_evolvability/WB053_3/N3_n=480_P.csv	Via & Conner	1995	WB053_3	N	P	480	2	180	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053_3.P	WB053_3.N
./data_evolvability/WB053b/A_n=80_G.csv	Via & Conner	1995	WB053b	A	G	80	2	182	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b.G	WB053b.A
./data_evolvability/WB053b/A_n=480_P.csv	Via & Conner	1995	WB053b	A	P	480	2	181	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b.P	WB053b.A
./data_evolvability/WB053b/N1_n=80_G.csv	Via & Conner	1995	WB053b	N	G	80	2	184	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b.G	WB053b.N
./data_evolvability/WB053b/N1_n=480_P.csv	Via & Conner	1995	WB053b	N	P	480	2	183	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b.P	WB053b.N
./data_evolvability/WB053b_3/A_n=80_G.csv	Via & Conner	1995	WB053b_3	A	G	80	2	186	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b_3.G	WB053b_3.A
./data_evolvability/WB053b_3/A_n=480_P.csv	Via & Conner	1995	WB053b_3	A	P	480	2	185	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b_3.P	WB053b_3.A
./data_evolvability/WB053b_3/N3_n=80_G.csv	Via & Conner	1995	WB053b_3	N	G	80	2	188	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b_3.G	WB053b_3.N
./data_evolvability/WB053b_3/N3_n=480_P.csv	Via & Conner	1995	WB053b_3	N	P	480	2	187	Tribolium	castaneum	half-sib	half-sib	insect	AN	WB053b_3.P	WB053b_3.N
./data_evolvability/WB056/A_n=64_G.csv	Service & Rose	1985	WB056	A	G	64	2	189	Drosophila	melanogaster	half-sib	half-sib	insect	AN	WB056.G	WB056.A
./data_evolvability/WB056/A_n=758_P.csv	Service & Rose	1985	WB056	A	P	758	2	190	Drosophila	melanogaster	half-sib	half-sib	insect	AN	WB056.P	WB056.A
./data_evolvability/WB056/N1_n=64_G.csv	Service & Rose	1985	WB056	N	G	64	2	191	Drosophila	melanogaster	half-sib	half-sib	insect	AN	WB056.G	WB056.N
./data_evolvability/WB056/N1_n=758_P.csv	Service & Rose	1985	WB056	N	P	758	2	192	Drosophila	melanogaster	half-sib	half-sib	insect	AN	WB056.P	WB056.N
./data_evolvability/WB057/A_n=22_G.csv	Grill et al.	1997	WB057	A	G	22	4	193	Harmonia	axyridis	full-sib	full-sib	insect	AN	WB057.G	WB057.A
./data_evolvability/WB057/A_n=277_P.csv	Grill et al.	1997	WB057	A	P	277	4	194	Harmonia	axyridis	full-sib	full-sib	insect	AN	WB057.P	WB057.A
./data_evolvability/WB057/N1_n=26_G.csv	Grill et al.	1997	WB057	N	G	26	4	195	Harmonia	axyridis	full-sib	full-sib	insect	AN	WB057.G	WB057.N
./data_evolvability/WB057/N1_n=277_P.csv	Grill et al.	1997	WB057	N	P	277	4	196	Harmonia	axyridis	full-sib	full-sib	insect	AN	WB057.P	WB057.N

We can then import the matrices and descriptive data.

         mats_evolve <- import_mat(matrices_evolve)
        dscrp_evolve <- import_mat(descripDat_evolve)

Scale the data to prevent disproportional large effects of the traits with large values:

  standard <- list()
  for (i in 1:length(dscrp_evolve)){
    # standardize descriptive data
    # calculate the average mean for each trait
    dscrp_evolve[[i]]$stanfac <- ave(dscrp_evolve[[i]]$mean,dscrp_evolve[[i]]$trait)
    #standardize mean, sd, se, Va
    dscrp_evolve[[i]]$mean <- dscrp_evolve[[i]]$mean/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$sd <- dscrp_evolve[[i]]$sd/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$se <- dscrp_evolve[[i]]$se/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$Va <- dscrp_evolve[[i]]$Va/(dscrp_evolve[[i]]$stanfac^2)
    # save the standardization factor in to another list
    standard[[i]] <- data.frame(trait=unique(dscrp_evolve[[i]]$trait),stanfac=dscrp_evolve[[i]]$stanfac[which(!duplicated(dscrp_evolve[[i]]$trait))]) 
    # remove spaces and dashes from the trait names to be consistent with the names in the covariance-matrices.
    standard[[i]]$trait <- gsub(" ","\\.", as.character(standard[[i]]$trait))
    standard[[i]]$trait <- gsub("-","\\.", standard[[i]]$trait)
    rownames(standard[[i]]) <- standard[[i]]$trait
    }
  
  id_d <- sapply(names(dscrp_evolve),function(x){strsplit(x,"/")[[1]][3]})
  id_m <- sapply(names(mats_evolve),function(x){strsplit(x,"/")[[1]][3]})
  for (i in 1:length(mats_evolve)){
    # select the right standardization factors
    j <- which(id_d == id_m[[i]])
    std <- standard[[j]]
    # remove traits not represented in the covar-matrix
    std <- std[which(std$trait %in% names(mats_evolve[[i]])),]
    mats_evolve[[i]] <- mats_evolve[[i]][which(names(mats_evolve[[i]]) %in% std$trait),]
    mats_evolve[[i]] <- mats_evolve[[i]][,which(names(mats_evolve[[i]]) %in% std$trait)]
    # order traits as in the covar-matrix
    std <- std[names(mats_evolve[[i]]),]
    # check whether the standardization factors line up
    if(sum(names(mats_evolve[[i]]) != std$trait)>0){stop("There is a inconsistency between trait names of a matrix and a description file")}
    mats_evolve[[i]] <- mats_evolve[[i]]/(std$stanfac %*% t(std$stanfac))
    }
  
  length(mats_evolve)

## [1] 196

  length(dscrp_evolve)

## [1] 49

2.2 – Descriptors of the Studies

Using these matrices we can also get some descriptor information about the studies

           env_mats_evolve <- extract_matrix_env(names(mats_evolve), type = "full")
         names_mats_evolve <- extract_matrix_studynames(names(mats_evolve))
          type_mats_evolve <- extract_matrix_type(names(mats_evolve))

        means_env <- plyr::ldply(lapply(dscrp_evolve, function(x) x[,c("env","trait","mean")]))
        
        # make trait names consistent with the row and column names
        means_env$trait <- gsub(" ","\\.",means_env$trait)
        means_env$trait <- gsub("_","\\.",means_env$trait)
        
    means_env$.id <- gsub("/DecsrpData.csv", "", means_env$.id)
    means_env$.id <- gsub("./data/", "", means_env$.id)
  means_env$env_r <- gsub("[1-9]", "", as.character(means_env$env))
  means_env$match <- paste0(means_env$.id, "_", means_env$env)

Now that we have extracted the mean vectors for all traits within each study and environment we can extract the G matrices for novel and non-novel environments.

           G_matrices <- mats_evolve[ grep("G", names(mats_evolve)) ]
        G_matrices_NS <- G_matrices[ grep("/[NS]", names(G_matrices)) ]
        G_matrices_AB <- G_matrices[ grep("/[AB]", names(G_matrices)) ]

        n_G_NS <- match_n(names(G_matrices_NS))
        n_G_AB <- match_n(names(G_matrices_AB))

We can check that everything matches, so we can be sure that nothing has been disorganized using some simple tests

    # First, we can subset just the G matrices from the metadata file
        tst <- meta_evolve %>% filter(type == "G")

    # Extract n from the G matrices
        n_tst <- match_n(names(G_matrices))

    # We can test that the sample sizes match what is actually in the file names
        all(tst$n == n_tst)  # TRUE, so we are good here.

## [1] FALSE

    # We can also make sure that the col matches with the G-matrices
        all(tst$col == match(tst$file, names(mats_evolve))) # TRUE, so good.

## [1] FALSE

    # Check that the sample sizes match with the matrices for each matrix. That way we can be absolutely sure that N matches the matrix
        n_NS_tst <- match_n(names(G_matrices_NS))
        n_AB_tst <- match_n(names(G_matrices_AB))

        all(n_G_NS == n_NS_tst) # TRUE

## [1] TRUE

        all(n_G_AB == n_AB_tst) # TRUE

## [1] TRUE

2.3 – Simulations for Evolvability Metrics

We now have 5000 simulated G matrices with the sample sizes (number of families) used to estimate these matrices. Given that we don’t have the P matrices for these data, we can’t simulate changes in the means using the same approaches (but see below). For now, we will use the mean differences in traits across environments to superimpose the plasticity vector on each G matrix. First, we need the centroid differences in means between the two environments:

    # To make things easy to compare, just split the data
    means_env$study <- gsub("./data_evolvability/", "", means_env$.id)
    split_mean_stdy <- base::split(means_env, means_env$study)
    split_mean_stdy <- lapply(split_mean_stdy, function(x) arrange(x, env))

    # Get the plasticity response vector between the non-novel / benign and novel / stressful.
    cen.diff <- lapply(split_mean_stdy, function(x)centroid.diff(x[x$env_r == "A" | x$env_r == "B",]$mean, x[x$env_r == "N" | x$env_r == "S",]$mean)$vector)
    
    # Get corresponding trait names
    cen.diff.names <- lapply(split_mean_stdy, function(x){as.character(x[x$env_r == "A" | x$env_r == "B",]$trait)})

Locate cases in which the number of traits do not match up

    ids <- unique(c(which(unlist(lapply(cen.diff,length))!=unlist(lapply(G_matrices_AB,function(x){dim(x)[1]}))),which(unlist(lapply(cen.diff,length))!=unlist(lapply(G_matrices_NS,function(x){dim(x)[1]})))))

    for(i in ids){
      # Identify the trait which all have in common
      intsect <- intersect(cen.diff.names[[i]],intersect(names(G_matrices_NS[[i]]),names(G_matrices_AB[[i]])))
      # select only the traits in common
      cen.diff[[i]] <- cen.diff[[i]][which(cen.diff.names[[i]]%in%intsect)]
      G_matrices_AB[[i]] <- G_matrices_AB[[i]][which(names(G_matrices_AB[[i]])%in%intsect),]
      G_matrices_AB[[i]] <- G_matrices_AB[[i]][,which(names(G_matrices_AB[[i]])%in%intsect)]
      G_matrices_NS[[i]] <- G_matrices_NS[[i]][which(names(G_matrices_NS[[i]])%in%intsect),]
      G_matrices_NS[[i]] <- G_matrices_NS[[i]][,which(names(G_matrices_NS[[i]])%in%intsect)]  
    }

Now that we have checked things and we have the novel (i.e., NS) and non-novel (i.e., AB) matrices extracted we can do some simulations with the matrices

    # Simulate G matrices for environments
        sims = 5000

        G_sim_NS <- mapply(function(x,y) sim_G(x, y, sims = sims, cor = FALSE, matrix_fix = "bend"), 
                            x = G_matrices_NS, 
                            y = split(n_G_NS, 1:length(n_G_NS)))
        G_sim_AB <- mapply(function(x,y) sim_G(x, y, sims = sims, cor = FALSE, matrix_fix = "bend"), 
                            x = G_matrices_AB, 
                            y = split(n_G_AB, 1:length(n_G_AB)))

Plasticity vectors are based on the multivariate mean centroid differences. We can now project these vectors on the G matrices in the non-novel and novel environment to understand how plastic responses brought about by the environment lead to greater, lower or no change in evolvability. We note that the trait means themselves also have sampling variances associated with them, however, we are assuming that these are fixed given we don’t have a P matrix to simulate the mean response vector for all studies. However, given that we are incorporating sampling error associated with G this will drive weighted differences among studies. Adding any additional sampling variance resulting from P vectors is likely not to change results because this will just add proportionally more sampling variance across studies in many cases.

    # Here we can get the amount of genetic variance along the plastic vector
        evolv_max_G_NS <- evolveability_along_plastic_vector(G_sim_NS, cen.diff, compare = "max")

        evolv_max_G_AB <- evolveability_along_plastic_vector(G_sim_AB, cen.diff, compare = "max")

    # Here we test whether the amount of genetic variance along the plastic vector is higher than expected
        evolv_max_G_NS_test <- evolveability_along_plastic_vector(G_sim_NS, cen.diff, compare = "maxcorr")

        evolv_max_G_AB_test <- evolveability_along_plastic_vector(G_sim_AB, cen.diff, compare = "maxcorr")
        
    # Here we can get the angle between the plastic vector and gmax
        evol_angle_NS <- evolveability_angle_along_plastic_vector(G_sim_NS, cen.diff)
        
        evol_angle_AB <- evolveability_angle_along_plastic_vector(G_sim_AB, cen.diff)

2.4 – Plotting Simulations

Now that we have all the simulations finished we can explore the distribution of effects across the studies

    plot_effect_dist(evolv_max_G_NS)

Figure 2.4.1 – Distribution of effects for the proportion of genetic variation along plastic response vector for \(G_{novel}\).

    plot_effect_dist(evolv_max_G_AB)

Figure 2.4.2 – Distribution of effects for the proportion of genetic variation along plastic response vector for \(G_{non-novel}\).

    plot_effect_dist(evol_angle_NS)

Figure 2.4.3 – Distribution of effects for the angle between \(gmax_{novel}\) and the plastic response vector for \(G_{novel}\).

    plot_effect_dist(evol_angle_AB)

Figure 2.4.4 – Distribution of effects for the angle between \(gmax_{nonnovel}\) and the plastic response vector for \(G_{non-novel}\).

2.5 – Generating the Data for Meta-analysis

We now want to condense down the simulated distributions and sampling variance so that we can more easily meta-analyse these data. Of course, one could technically use the entire distribution, but computationally this will make the models take a long time to run. So, we will keep it simple. We will first generate the data for \(\pi_{e}\).

                effect_G_NS <- create_data_evolve(evolv_max_G_NS, meta_evolve, env = c("N", "S"))

                effect_G_AB <- create_data_evolve(evolv_max_G_AB, meta_evolve, env = c("A", "B"))

            evolve_max_data <- rbind(effect_G_NS, effect_G_AB)
        evolve_max_data$err <- 1:nrow(evolve_max_data)
evolve_max_data$environment <- ifelse(evolve_max_data$env == "B" | evolve_max_data$env == "A", "Non-novel", "Novel")
         evolve_max_data$stdy <- substring(evolve_max_data$study,1, 5)
       
        write.csv(evolve_max_data, file = "./data_evolvability/evolve_max_data_standardised.csv",row.names=FALSE)

Next for will generate the data for \(\pi_{e-0}\).

                effect_G_NS_test <- create_data_evolve(evolv_max_G_NS_test, meta_evolve, env = c("N", "S"))

                effect_G_AB_test <- create_data_evolve(evolv_max_G_AB_test, meta_evolve, env = c("A", "B"))

            evolve_max_data_test <- rbind(effect_G_NS_test, effect_G_AB_test)
        evolve_max_data_test$err <- 1:nrow(evolve_max_data_test)
evolve_max_data_test$environment <- ifelse(evolve_max_data_test$env == "B" | evolve_max_data_test$env == "A", "Non-novel", "Novel")
         evolve_max_data_test$stdy <- substring(evolve_max_data_test$study,1, 5)
       
        write.csv(evolve_max_data_test, file = "./data_evolvability/evolve_max_data_test.csv",row.names=FALSE)

Now that we have data for \(pi_{e}\), we can generate \(\theta_{e}\) in the same way:

                     angle_G_NS <- create_data_evolve(evol_angle_NS, meta_evolve, env = c("N", "S"))

                     angle_G_AB <- create_data_evolve(evol_angle_AB, meta_evolve, env = c("A", "B"))

              evolve_angle_data <- rbind(angle_G_NS, angle_G_AB)
          evolve_angle_data$err <- 1:nrow(evolve_angle_data)
  evolve_angle_data$environment <- ifelse(evolve_angle_data$env == "B" | evolve_angle_data$env == "A", "Non-novel", "Novel")
       evolve_angle_data$stdy <- substring(evolve_angle_data$study,1, 5)


        write.csv(evolve_angle_data, file = "./data_evolvability/evolve_angle_data_standardised.csv", row.names=FALSE)

3.0 – Meta-analysis 1: General patterns of G and P within and across environments

3.1 – Load Data

    # Load libraries
        library(metafor)
        library(knitr)
        library(rmarkdown)

    # Load data
                SMD_data <- read.csv("./data_GP_compare/SMD_data_standardised.csv")
              angle_data <- read.csv("./data_GP_compare/angle_across_data_standardised.csv")
       angle_data_within <- read.csv("./data_GP_compare/angle_within_data_standardised.csv")
    var_along_logit_data <- read.csv("./data_GP_compare/var_along_logit_data_standardised.csv")

3.2 – Meta-analysis of SDV

First, we can have a look at the dataset more carefully before running the analysis.

    kable(SMD_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

X	study	type	ES	S_var	authors	year	avg_n_type	trait_num	genus	species	design	group	env.comp	species_full	stdy	err	wi
1	RR073	G	-1.1114684	0.0137012	Evans et al.	2015	66.0	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	1	72.98608
2	RR073	P	-1.1761978	0.0036433	Evans et al.	2015	153.5	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	2	274.47870
3	RR076	G	-1.3864048	0.0141709	Guan et al.	2016	69.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	3	70.56698
4	RR076	P	-0.4165704	0.0010826	Guan et al.	2016	2525.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	4	923.73149
5	RR087a	G	-0.0463060	0.0266880	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	5	37.47005
6	RR087a	P	-0.1584050	0.0011757	Sae-Lim et al.	2013	2127.5	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	6	850.52360
7	RR087b	G	-0.0811399	0.0251500	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	7	39.76148
8	RR087b	P	0.0924924	0.0009454	Sae-Lim et al.	2013	2577.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	8	1057.73026
9	RR087c	G	-0.2067433	0.0217486	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	9	45.97998
10	RR087c	P	-0.0052831	0.0010939	Sae-Lim et al.	2013	2091.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	10	914.14934
11	sR097	G	-1.0757071	0.0647484	Kasule	1991	26.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	11	15.44440
12	sR097	P	0.5490437	0.0073155	Kasule	1991	312.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	12	136.69640
13	sR110b	G	0.2457675	0.0971304	Blanckenhorn & Heyland	2004	32.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	13	10.29544
14	sR110b	P	0.1302658	0.0254568	Blanckenhorn & Heyland	2004	113.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	14	39.28231
15	sR110c	G	-0.1644129	0.0885876	Blanckenhorn & Heyland	2004	36.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	15	11.28826
16	sR110c	P	-0.2031563	0.0200968	Blanckenhorn & Heyland	2004	150.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	16	49.75914
17	sR115	G	-0.9262440	0.0020576	Conner et al.	2003	471.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	17	485.99435
18	sR115	P	0.4498877	0.0017084	Conner et al.	2003	1010.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	18	585.34934
19	sR116	G	-0.3003715	0.0069157	Delcourt et al.	2010	365.5	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	19	144.59791
20	sR116	P	-0.0686755	0.0009931	Delcourt et al.	2010	3299.0	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	20	1006.93838
21	sR145	G	0.1496778	0.0308049	Tucic et al.	1991	47.5	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	21	32.46242
22	sR145	P	0.0567341	0.0116437	Tucic et al.	1991	137.0	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	22	85.88331
23	WB002a	G	-0.0332739	0.0675742	Begin & Roff	2001	39.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	23	14.79854
24	WB002a	P	-0.2989204	0.0033918	Begin & Roff	2001	605.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	24	294.82803
25	WB003	G	0.8505977	0.0343841	Begin et al.	2004	59.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	25	29.08320
26	WB003	P	0.4246671	0.0038752	Begin et al.	2004	697.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	26	258.05423
27	WB003_2	G	0.3356583	0.0435435	Begin et al.	2004	58.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	27	22.96552
28	WB003_2	P	0.3286113	0.0043748	Begin et al.	2004	624.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	28	228.58439
29	WB003a	G	-0.0378720	0.0682350	Begin et al.	2004	43.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	29	14.65523
30	WB003a	P	0.3435527	0.0093449	Begin et al.	2004	288.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	30	107.01026
31	WB003a_2	G	0.7981247	0.0636134	Begin et al.	2004	38.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	31	15.71997
32	WB003a_2	P	0.4089151	0.0164377	Begin et al.	2004	197.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	32	60.83565
33	WB007	G	-1.0716421	0.0068572	Brock et al.	2010	131.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	33	145.83305
34	WB007	P	-1.1222659	0.0011237	Brock et al.	2010	562.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	34	889.90631
35	WB010	G	0.1409895	0.0374687	Collins et al.	1999	101.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	35	26.68894
36	WB010	P	0.0187338	0.0042686	Collins et al.	1999	991.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	36	234.26914
37	WB011	G	0.2529621	0.0095757	Czesak & Fox	2003	404.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	37	104.43115
38	WB011	P	0.0000158	0.0010038	Czesak & Fox	2003	2667.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	38	996.20709
39	WB028	G	1.3743055	0.0162181	Klause & Morin	2001	60.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	39	61.65935
40	WB028	P	1.4176906	0.0012685	Klause & Morin	2001	600.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	40	788.30381
41	WB029	G	-0.2907661	0.0601083	King et al.	2011	63.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	41	16.63664
42	WB029	P	-1.0680634	0.0008849	King et al.	2011	1695.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	42	1130.05446
43	WB033	G	0.0528833	0.0663470	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	43	15.07227
44	WB033	P	0.5024050	0.0069455	Lau et al.	2014	513.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	44	143.97718
45	WB033a	G	0.9189091	0.0414328	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	45	24.13545
46	WB033a	P	-0.0210253	0.0089497	Lau et al.	2014	464.5	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	46	111.73533
47	WB035	G	-0.0573465	0.0064297	Messina et al.	2013	432.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	47	155.52774
48	WB035	P	-1.2368324	0.0004308	Messina et al.	2013	2204.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	48	2321.03330
49	WB038	G	-0.0589241	0.0177160	Rauter&Moore	2002	135.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	49	56.44611
50	WB038	P	-0.3677845	0.0032270	Rauter&Moore	2002	639.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	50	309.88875
51	WB045	G	1.6727335	0.0085023	Simons & Roff	1996	39.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	51	117.61495
52	WB045	P	0.2529547	0.0076488	Simons & Roff	1996	371.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	52	130.73914
53	WB045a	G	0.5721094	0.0636117	Simons & Roff	1996	39.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	53	15.72038
54	WB045a	P	0.1596246	0.0067022	Simons & Roff	1996	381.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	54	149.20430
55	WB053	G	-1.3688710	0.0476525	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	55	20.98524
56	WB053	P	-0.0174303	0.0055034	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	56	181.70564
57	WB053_3	G	-1.3669331	0.0451303	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	57	22.15807
58	WB053_3	P	0.0035190	0.0048829	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	58	204.79462
59	WB053b	G	-1.3857924	0.0135062	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	59	74.04016
60	WB053b	P	-0.0149641	0.0060127	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	60	166.31377
61	WB053b_3	G	-1.3855821	0.0116041	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	61	86.17632
62	WB053b_3	P	0.0045548	0.0064249	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	62	155.64431
63	WB057	G	1.1158349	0.0766307	Grill et al.	1997	24.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	63	13.04959
64	WB057	P	0.4308134	0.0117949	Grill et al.	1997	277.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	64	84.78246

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual and between-study variance. We can also can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Meta-analysis - Intercept Only Model
        # G 
            G_SMD_data <- subset(SMD_data, type == "G")
            
            model_intcp_SMD_het_G <- rma.mv(ES ~ 1, V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)

        # P 
            P_SMD_data <- subset(SMD_data, type == "P")
            
            model_intcp_SMD_het_P <- rma.mv(ES ~ 1, V = P_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = P_SMD_data)

    # Extract the variance estimates from the model
        #G
             stdy_sig_SMD_G <- model_intcp_SMD_het_G$sigma2[1]
              err_sig_SMD_G <- model_intcp_SMD_het_G$sigma2[2]

        #P
             stdy_sig_SMD_P <- model_intcp_SMD_het_P$sigma2[1]
              err_sig_SMD_P <- model_intcp_SMD_het_P$sigma2[2]

    # Calculate total sampling variance
        #G
                wi_SMD_G <- 1 / G_SMD_data$S_var
                vi_SMD_G <- sum(wi_SMD_G) * ( nrow(G_SMD_data) - 1) / ( (sum(wi_SMD_G)^2) - sum(wi_SMD_G))
                
        #P
                wi_SMD_P <- 1 / P_SMD_data$S_var
                vi_SMD_P <- sum(wi_SMD_P) * ( nrow(P_SMD_data) - 1) / ( (sum(wi_SMD_P)^2) - sum(wi_SMD_P))

    # Generate points estimates for I2 – heterogeneity at different levels
        #G
            I2total_SMD_G <- round((stdy_sig_SMD_G + err_sig_SMD_G) / (stdy_sig_SMD_G + err_sig_SMD_G + vi_SMD_G), digits = 2)

            I2stdy_SMD_G <- round((stdy_sig_SMD_G) / (stdy_sig_SMD_G + err_sig_SMD_G + vi_SMD_G), digits = 2)

        #P
            I2total_SMD_P <- round((stdy_sig_SMD_P + err_sig_SMD_P) / (stdy_sig_SMD_P + err_sig_SMD_P + vi_SMD_P), digits = 2)

            I2stdy_SMD_P <- round((stdy_sig_SMD_P) / (stdy_sig_SMD_P + err_sig_SMD_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_SMD = 0.98 and between study heterogeneity as I2stdy_SMD = 0.84) and for P (I2total_SMD = 0.99 and between study heterogeneity as I2stdy_SMD = 0.97). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should use a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actual effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this will reflect the median of the distribution of untransformed values rather than the mean.

    # Meta-regression models
        model_SMD_G <- metafor::rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        model_SMD_G_stress <- metafor::rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        model_SMD_G_half <- metafor::rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        
         n_design_SMD <- table(G_SMD_data$design)
        n_encComp_smd <- table(G_SMD_data$env.comp)
        
        # Marginalised means for different groups
        SMD_D_AN <- marginal(model_SMD_G, "design", G_SMD_data)
        SMD_D_BS <- marginal(model_SMD_G_stress, "design", G_SMD_data)

        AN_BS_SMD <- cbind(rbind(SMD_D_AN, SMD_D_BS), n_encComp_smd)

        SMD_D_full <- marginal(model_SMD_G, "env.comp", G_SMD_data)
        SMD_D_half <- marginal(model_SMD_G_half, "env.comp", G_SMD_data)

        design_SMD <- cbind(rbind(SMD_D_full, SMD_D_half), n_design_SMD)

        SMD_marg_table <- rbind(AN_BS_SMD, design_SMD)

        # We will just fit the P model
        model_SMD_P <- metafor::rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = P_SMD_data)

Once we have these models run we can have a look at the table of estimates:

    table_SMD <- data.frame(Est. = round(rbind(model_SMD_G$b, model_SMD_P$b), digits =3), L95CI = round(c(model_SMD_G$ci.lb, model_SMD_P$ci.lb), digits=3), U95CI = round(c(model_SMD_G$ci.ub, model_SMD_P$ci.ub), digits =3))
    rownames(table_SMD) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_SMD) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept_G	0.623	0.086	1.160
Half-sib Breeding Design_G	-1.092	-1.751	-0.434
Stressful Environment_G	-0.011	-0.633	0.612
Trait Number_G	-0.090	-0.359	0.179
Intercept_P	0.444	0.011	0.877
Half-sib Breeding Design_P	-0.531	-1.059	-0.002
Stressful Environment_P	-0.310	-0.809	0.190
Trait Number_P	-0.200	-0.411	0.011

3.3 – Meta-analysis of Proportion of Variance Along max

First, we can have a look at the dataset more carefully before running the analysis.

    kable(var_along_logit_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

X	study	type	ES	S_var	authors	year	avg_n_type	trait_num	genus	species	design	group	env.comp	species_full	stdy	err	wi
1	RR073	G	0.3854686	0.0881714	Evans et al.	2015	66.0	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	1	11.341545
2	RR073	P	0.4469213	0.0303974	Evans et al.	2015	153.5	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	2	32.897531
3	RR076	G	-0.9831790	0.1172049	Guan et al.	2016	69.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	3	8.532065
4	RR076	P	-0.3773049	0.0031580	Guan et al.	2016	2525.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	4	316.654156
5	RR087a	G	0.9952143	0.0739473	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	5	13.523150
6	RR087a	P	0.2297012	0.0037580	Sae-Lim et al.	2013	2127.5	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	6	266.100828
7	RR087b	G	0.4581764	0.0803353	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	7	12.447823
8	RR087b	P	-0.0555089	0.0031329	Sae-Lim et al.	2013	2577.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	8	319.193821
9	RR087c	G	-0.0984979	0.0750109	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	9	13.331403
10	RR087c	P	-0.5183770	0.0038306	Sae-Lim et al.	2013	2091.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	10	261.055378
11	sR097	G	-1.1636097	0.3387223	Kasule	1991	26.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	11	2.952271
12	sR097	P	-0.0071598	0.0253223	Kasule	1991	312.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	12	39.490913
13	sR110b	G	-2.8170762	0.2759173	Blanckenhorn & Heyland	2004	32.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	13	3.624274
14	sR110b	P	-0.5508976	0.0727283	Blanckenhorn & Heyland	2004	113.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	14	13.749807
15	sR110c	G	-1.0841059	0.2345270	Blanckenhorn & Heyland	2004	36.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	15	4.263901
16	sR110c	P	-0.6293073	0.0537564	Blanckenhorn & Heyland	2004	150.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	16	18.602440
17	sR115	G	0.3037840	0.0127205	Conner et al.	2003	471.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	17	78.613406
18	sR115	P	1.0542591	0.0051981	Conner et al.	2003	1010.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	18	192.378125
19	sR116	G	-0.3731776	0.0201041	Delcourt et al.	2010	365.5	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	19	49.741126
20	sR116	P	-0.0689952	0.0016041	Delcourt et al.	2010	3299.0	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	20	623.388723
21	sR145	G	0.5545398	0.0964286	Tucic et al.	1991	47.5	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	21	10.370367
22	sR145	P	-0.1463109	0.0374098	Tucic et al.	1991	137.0	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	22	26.730963
23	WB002a	G	0.0857084	0.1747969	Begin & Roff	2001	39.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	23	5.720927
24	WB002a	P	-0.2802787	0.0100423	Begin & Roff	2001	605.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	24	99.578399
25	WB003	G	0.5040949	0.0973839	Begin et al.	2004	59.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	25	10.268643
26	WB003	P	0.4346733	0.0086351	Begin et al.	2004	697.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	26	115.806945
27	WB003_2	G	0.0018268	0.0991137	Begin et al.	2004	58.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	27	10.089421
28	WB003_2	P	0.2619965	0.0090460	Begin et al.	2004	624.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	28	110.545546
29	WB003a	G	0.2298057	0.1362611	Begin et al.	2004	43.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	29	7.338851
30	WB003a	P	0.6737641	0.0191354	Begin et al.	2004	288.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	30	52.259231
31	WB003a_2	G	1.0096968	0.1659049	Begin et al.	2004	38.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	31	6.027551
32	WB003a_2	P	0.8927405	0.0305538	Begin et al.	2004	197.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	32	32.729101
33	WB007	G	0.7948070	0.0389207	Brock et al.	2010	131.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	33	25.693256
34	WB007	P	0.0222790	0.0089746	Brock et al.	2010	562.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	34	111.425492
35	WB010	G	-0.1958279	0.0877272	Collins et al.	1999	101.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	35	11.398971
36	WB010	P	-0.2305906	0.0103322	Collins et al.	1999	991.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	36	96.784861
37	WB011	G	1.4972372	0.0200357	Czesak & Fox	2003	404.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	37	49.910840
38	WB011	P	0.2038146	0.0029356	Czesak & Fox	2003	2667.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	38	340.646298
39	WB028	G	2.4763329	0.1339905	Klause & Morin	2001	60.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	39	7.463214
40	WB028	P	2.4941067	0.0130883	Klause & Morin	2001	600.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	40	76.404187
41	WB029	G	0.1310828	0.1289894	King et al.	2011	63.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	41	7.752577
42	WB029	P	0.1849028	0.0041783	King et al.	2011	1695.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	42	239.332622
43	WB033	G	-0.1928236	0.1202423	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	43	8.316539
44	WB033	P	0.5169502	0.0128358	Lau et al.	2014	513.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	44	77.907073
45	WB033a	G	0.7918906	0.1318105	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	45	7.586650
46	WB033a	P	0.0045761	0.0141276	Lau et al.	2014	464.5	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	46	70.783494
47	WB035	G	-1.4463889	0.0149326	Messina et al.	2013	432.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	47	66.967617
48	WB035	P	-2.5735619	0.0030416	Messina et al.	2013	2204.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	48	328.770705
49	WB038	G	0.2309286	0.0439186	Rauter&Moore	2002	135.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	49	22.769383
50	WB038	P	-0.5752013	0.0091143	Rauter&Moore	2002	639.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	50	109.717690
51	WB045	G	2.8750053	0.1531809	Simons & Roff	1996	39.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	51	6.528228
52	WB045	P	0.5532212	0.0149558	Simons & Roff	1996	371.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	52	66.863602
53	WB045a	G	0.5798272	0.1588199	Simons & Roff	1996	39.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	53	6.296442
54	WB045a	P	0.2880374	0.0152644	Simons & Roff	1996	381.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	54	65.511946
55	WB053	G	-2.7022412	0.3643009	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	55	2.744983
56	WB053	P	0.2136209	0.0165484	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	56	60.428861
57	WB053_3	G	-4.7140006	0.3802670	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	57	2.629731
58	WB053_3	P	-0.2077135	0.0165950	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	58	60.259201
59	WB053b	G	-2.8298645	0.1024369	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	59	9.762108
60	WB053b	P	-1.0942523	0.0163428	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	60	61.189195
61	WB053b_3	G	-4.7732524	0.1007585	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	61	9.924717
62	WB053b_3	P	-0.1414758	0.0169504	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	62	58.995497
63	WB057	G	-0.4329326	0.3313705	Grill et al.	1997	24.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	63	3.017770
64	WB057	P	-0.5492158	0.0249628	Grill et al.	1997	277.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	64	40.059601

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity of effects. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Multi-level meta-analysis
    # Subset the data and run the intercept only models
        # G
            G_varAlong_data <- subset(var_along_logit_data, type == "G")
            
            model_intcp_varAlong_het_G <- rma.mv(ES ~ 1, V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)

        # P
            P_varAlong_data <- subset(var_along_logit_data, type == "P")
            
            model_intcp_varAlong_het_P <- rma.mv(ES ~ 1, V = P_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = P_varAlong_data)

    # Extract the variance estimates from the model
        #G
            stdy_sig_varAlong_G <- model_intcp_varAlong_het_G$sigma2[1]
             err_sig_varAlong_G <- model_intcp_varAlong_het_G$sigma2[2]

        #P
             stdy_sig_varAlong_P <- model_intcp_varAlong_het_P$sigma2[1]
             err_sig_varAlong_P <- model_intcp_varAlong_het_P$sigma2[2]

    # Calculate total sampling variance
        #G
                wi_varAlong_G <- 1 / G_varAlong_data$S_var
                vi_varAlong_G <- sum(wi_varAlong_G) * ( nrow(G_varAlong_data) - 1) / ( (sum(wi_varAlong_G)^2) - sum(wi_varAlong_G))
        #P
                wi_varAlong_P <- 1 / P_varAlong_data$S_var
                vi_varAlong_P <- sum(wi_varAlong_P) * ( nrow(P_varAlong_data) - 1) / ( (sum(wi_varAlong_P)^2) - sum(wi_varAlong_P))

    # Generate points estimates for I2 – heterogenity at different levels
        #G
            I2total_varAlong_G <- round((stdy_sig_varAlong_G + err_sig_varAlong_G) / (stdy_sig_varAlong_G + err_sig_varAlong_G + vi_varAlong_G), digits = 2)

            I2stdy_varAlong_G <- round((stdy_sig_varAlong_G) / (stdy_sig_varAlong_G + err_sig_varAlong_G + vi_SMD_G), digits = 2)

        # P
            I2total_varAlong_P <- round((stdy_sig_varAlong_P + err_sig_varAlong_P) / (stdy_sig_varAlong_P + err_sig_varAlong_P + vi_varAlong_P), digits = 2)

            I2stdy_varAlong_P <- round((stdy_sig_varAlong_P) / (stdy_sig_varAlong_P + err_sig_varAlong_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_Varlong = 0.97 and between study heterogeneity as I2stdy_varAlong = 0.73) and for P (I2total_varAlong = 0.99 and between study heterogeneity as I2stdy_varAlong = 0.82). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

   # Total variance change along g/p max.
        model_VarAlong_G <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)
        model_VarAlong_G_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)
        model_VarAlong_G_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)

        n_AN_S <- table(G_varAlong_data$env.comp)
        n_design <- table(G_varAlong_data$design)

        # Marginal means
        varalong_G_AN <- marginal(model_VarAlong_G, "design", data = G_varAlong_data)
        varalong_G_S  <- marginal(model_VarAlong_G_stress, "design", data = G_varAlong_data)

        env_cmp_varAlong <- rbind(varalong_G_AN, varalong_G_S)
        env_cmp_varAlong <- cbind(env_cmp_varAlong, n_AN_S)

        varalong_G_full <- marginal(model_VarAlong_G, "env.comp", G_varAlong_data)
        varalong_G_half <- marginal(model_VarAlong_G_half, "env.comp", G_varAlong_data)
          
        design_varAlong <- rbind(varalong_G_full, varalong_G_half)
        design_varAlong <- cbind(design_varAlong, n_design)

        varalong_marg_means <- rbind(env_cmp_varAlong, design_varAlong)

        model_VarAlong_P <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = P_varAlong_data)

Once we have these models run we can have a look at the table of estimates:

    table_varAlong <- data.frame(Est. = round(rbind(model_VarAlong_G$b, model_VarAlong_P$b), digits =3), L95CI = round(c(model_VarAlong_G$ci.lb, model_VarAlong_P$ci.lb), digits=3), U95CI = round(c(model_VarAlong_G$ci.ub, model_VarAlong_P$ci.ub), digits =3))
    rownames(table_varAlong) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_varAlong) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept_G	0.466	-0.501	1.432
Half-sib Breeding Design_G	-1.479	-2.681	-0.278
Stressful Environment_G	0.753	-0.387	1.893
Trait Number_G	0.453	-0.045	0.952
Intercept_P	0.320	-0.392	1.031
Half-sib Breeding Design_P	-0.517	-1.393	0.360
Stressful Environment_P	0.011	-0.819	0.841
Trait Number_P	0.070	-0.289	0.430

3.4 – Meta-analysis of Angle Differences Across Environments

First, we can have a look at the dataset more carefully before running the analysis.

    kable(angle_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

X	study	type	ES	S_var	authors	year	avg_n_type	trait_num	genus	species	design	group	env.comp	species_full	stdy	err	wi
1	RR073	G	0.2525296	0.6071997	Evans et al.	2015	66.0	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	1	1.6469047
2	RR073	P	1.9626238	0.4822250	Evans et al.	2015	153.5	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	2	2.0737207
3	RR076	G	-2.1552714	0.0548510	Guan et al.	2016	69.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	3	18.2311949
4	RR076	P	-3.0958212	0.0199065	Guan et al.	2016	2525.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	4	50.2348513
5	RR087a	G	-2.4982323	0.3015269	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	5	3.3164541
6	RR087a	P	-1.8113762	0.0168979	Sae-Lim et al.	2013	2127.5	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	6	59.1788134
7	RR087b	G	-2.4223537	0.7267958	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	7	1.3759022
8	RR087b	P	-2.2142142	0.0399846	Sae-Lim et al.	2013	2577.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	8	25.0096430
9	RR087c	G	-2.4953079	0.7500006	Sae-Lim et al.	2013	100.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	9	1.3333322
10	RR087c	P	-1.4461472	0.0402833	Sae-Lim et al.	2013	2091.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	10	24.8242095
11	sR097	G	-1.9352107	0.8905764	Kasule	1991	26.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	11	1.1228683
12	sR097	P	-2.7588449	0.6487915	Kasule	1991	312.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	12	1.5413271
13	sR110b	G	-1.3485285	0.4753950	Blanckenhorn & Heyland	2004	32.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	13	2.1035140
14	sR110b	P	-1.5643303	0.1360035	Blanckenhorn & Heyland	2004	113.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	14	7.3527511
15	sR110c	G	-2.6295399	1.2333280	Blanckenhorn & Heyland	2004	36.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	15	0.8108143
16	sR110c	P	-2.2186783	0.2409692	Blanckenhorn & Heyland	2004	150.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	16	4.1499077
17	sR115	G	-0.8877191	0.0241688	Conner et al.	2003	471.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	17	41.3756926
18	sR115	P	-1.2472474	0.0135478	Conner et al.	2003	1010.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	18	73.8126450
19	sR116	G	-3.4919763	0.2700806	Delcourt et al.	2010	365.5	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	19	3.7025989
20	sR116	P	-4.8715979	0.1116487	Delcourt et al.	2010	3299.0	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	20	8.9566685
21	sR145	G	1.5171679	1.2131713	Tucic et al.	1991	47.5	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	21	0.8242859
22	sR145	P	-2.2373713	0.1405605	Tucic et al.	1991	137.0	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	22	7.1143720
23	WB002a	G	-2.5297233	0.4065465	Begin & Roff	2001	39.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	23	2.4597430
24	WB002a	P	-2.9790205	0.1099553	Begin & Roff	2001	605.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	24	9.0946022
25	WB003	G	-2.6296686	0.1720626	Begin et al.	2004	59.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	25	5.8118366
26	WB003	P	-3.3843734	0.0551579	Begin et al.	2004	697.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	26	18.1297820
27	WB003_2	G	-2.5908291	0.1617689	Begin et al.	2004	58.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	27	6.1816566
28	WB003_2	P	-3.3986432	0.1291171	Begin et al.	2004	624.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	28	7.7449093
29	WB003a	G	-1.8234313	0.0993599	Begin et al.	2004	43.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	29	10.0644176
30	WB003a	P	-3.1296428	0.1030189	Begin et al.	2004	288.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	30	9.7069562
31	WB003a_2	G	-2.1081164	0.0584278	Begin et al.	2004	38.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	31	17.1151325
32	WB003a_2	P	-2.9168284	0.0972165	Begin et al.	2004	197.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	32	10.2863196
33	WB007	G	0.8783017	0.0321940	Brock et al.	2010	131.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	33	31.0616814
34	WB007	P	-1.4142582	0.0086498	Brock et al.	2010	562.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	34	115.6091924
35	WB010	G	-1.5891018	0.0166659	Collins et al.	1999	101.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	35	60.0028714
36	WB010	P	-3.5814923	0.3460956	Collins et al.	1999	991.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	36	2.8893750
37	WB011	G	-5.3388824	0.2657256	Czesak & Fox	2003	404.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	37	3.7632799
38	WB011	P	-2.2309386	0.0172158	Czesak & Fox	2003	2667.0	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	38	58.0862199
39	WB028	G	-0.2645351	0.0367901	Klause & Morin	2001	60.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	39	27.1812343
40	WB028	P	-0.6601235	0.0217321	Klause & Morin	2001	600.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	40	46.0148668
41	WB029	G	-4.3993865	0.6267768	King et al.	2011	63.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	41	1.5954643
42	WB029	P	-0.8048412	0.0024736	King et al.	2011	1695.0	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	42	404.2631712
43	WB033	G	-5.5539566	0.5108537	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	43	1.9575077
44	WB033	P	-6.4522773	0.3105651	Lau et al.	2014	513.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	44	3.2199364
45	WB033a	G	-4.6843816	0.0960867	Lau et al.	2014	60.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	45	10.4072693
46	WB033a	P	-6.1899906	0.1131624	Lau et al.	2014	464.5	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	46	8.8368550
47	WB035	G	-0.6864944	0.0175390	Messina et al.	2013	432.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	47	57.0157640
48	WB035	P	-1.7223657	0.0215781	Messina et al.	2013	2204.0	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	48	46.3432925
49	WB038	G	-2.0606936	0.2091481	Rauter&Moore	2002	135.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	49	4.7813013
50	WB038	P	-2.7114595	0.0360297	Rauter&Moore	2002	639.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	50	27.7548998
51	WB045	G	-0.2207482	0.0074963	Simons & Roff	1996	39.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	51	133.3985422
52	WB045	P	-3.0747363	0.1228861	Simons & Roff	1996	371.0	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	52	8.1376165
53	WB045a	G	-2.8350612	0.2934883	Simons & Roff	1996	39.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	53	3.4072906
54	WB045a	P	-2.9119533	0.1649474	Simons & Roff	1996	381.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	54	6.0625388
55	WB053	G	-2.8334129	0.4310894	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	55	2.3197045
56	WB053	P	-3.5804857	1.1341314	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	56	0.8817320
57	WB053_3	G	-2.5672827	1.1912009	Via & Conner	1995	24.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	57	0.8394889
58	WB053_3	P	-0.4454588	0.0396640	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	58	25.2117556
59	WB053b	G	-3.0341318	0.1633039	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	59	6.1235516
60	WB053b	P	-1.6949102	0.0586645	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	60	17.0460785
61	WB053b_3	G	-3.1392771	1.1214314	Via & Conner	1995	80.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	61	0.8917175
62	WB053b_3	P	-3.0306846	0.2700453	Via & Conner	1995	480.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	62	3.7030825
63	WB057	G	-1.9043333	0.1313720	Grill et al.	1997	24.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	63	7.6119738
64	WB057	P	-3.1128513	0.0168762	Grill et al.	1997	277.0	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	64	59.2549553

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Subset the data and run the intercept only models
     # G
        G_angle_data <- subset(angle_data, type == "G")
            
        model_intcp_angle_het_G <- rma.mv(ES ~ 1, V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)

    # P
        P_angle_data <- subset(angle_data, type == "P")
            
        model_intcp_angle_het_P <- rma.mv(ES ~ 1, V = P_angle_data$S_var, random = list(~1|stdy, ~1|err), data = P_angle_data)
            
    # Extract the variance estimates from the model     
        # G 
            stdy_sig_angle_G <- model_intcp_angle_het_G$sigma2[1]
             err_sig_angle_G <- model_intcp_angle_het_G$sigma2[2]
        # P
             stdy_sig_angle_P <- model_intcp_angle_het_P$sigma2[1]
             err_sig_angle_P <- model_intcp_angle_het_P$sigma2[2]

    # Calculate total sampling variance
        # G
                wi_angle_G <- 1 / G_angle_data$S_var
                vi_angle_G <- sum(wi_angle_G) * ( nrow(G_angle_data) - 1) / ( (sum(wi_angle_G)^2) - sum(wi_angle_G))
        # P 
                wi_angle_P <- 1 / P_angle_data$S_var
                vi_angle_P <- sum(wi_angle_P) * ( nrow(P_angle_data) - 1) / ( (sum(wi_angle_P)^2) - sum(wi_angle_P))
    
    # Generate points estimates for I2 – heterogenity at different levels
        #G
            I2total_angle_G <- round((stdy_sig_angle_G + err_sig_angle_G) / (stdy_sig_angle_G + err_sig_angle_G + vi_angle_G), digits = 2)

            I2stdy_angle_G <- round((stdy_sig_angle_G) / (stdy_sig_angle_G + err_sig_angle_G + vi_SMD_G), digits = 2)

        #P
            I2total_angle_P <- round((stdy_sig_angle_P + err_sig_angle_P) / (stdy_sig_angle_P + err_sig_angle_P + vi_angle_P), digits = 2)

            I2stdy_angle_P <- round((stdy_sig_angle_P) / (stdy_sig_angle_P + err_sig_angle_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_angle = 0.98 and between study heterogeneity as I2stdy_angle = 0.85) and for P (I2total_angle = 0.99 and between study heterogeneity as I2stdy_angle = 0.88). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

    # Multi-level meta-regression models of the angle changes for G and P across environments. 
               model_Angle_G <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)
        model_Angle_G_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)
          model_Angle_G_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)

        n_AN_S_angle <- table(G_angle_data$env.comp)
        n_design_angle <- table(G_angle_data$design)

        # Marginal means
        angle_G_AN <- marginal(model_Angle_G, "design", data = G_angle_data)
        angle_G_S  <- marginal(model_Angle_G_stress, "design", data = G_angle_data)

        env_cmp_angle <- rbind(angle_G_AN, angle_G_S)
        env_cmp_angle <- cbind(env_cmp_angle, n_AN_S_angle)

        Angle_G_full <- marginal(model_Angle_G, "env.comp", G_angle_data)
        Angle_G_half <- marginal(model_Angle_G_half, "env.comp", G_angle_data)
          
        design_angle <- rbind(Angle_G_full, Angle_G_half)
        design_angle <- cbind(design_angle, n_design_angle)

        angle_marg_means <- rbind(env_cmp_angle, design_angle)

          model_Angle_P <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_angle_data$S_var, random = list(~1|stdy, ~1|err), data = P_angle_data)

Once we have these models run we can have a look at the table of estimates:

    table_angle <- data.frame(Est. = round(rbind(model_Angle_G$b, model_Angle_P$b), digits =3), L95CI = round(c(model_Angle_G$ci.lb, model_Angle_P$ci.lb), digits=3), U95CI = round(c(model_Angle_G$ci.ub, model_Angle_P$ci.ub), digits =3))
    rownames(table_angle) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_angle) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept_G	-1.972	-3.015	-0.928
Half-sib Breeding Design_G	-1.246	-2.551	0.058
Stressful Environment_G	1.513	0.273	2.753
Trait Number_G	0.891	0.353	1.428
Intercept_P	-2.923	-4.062	-1.784
Half-sib Breeding Design_P	-0.511	-1.920	0.898
Stressful Environment_P	1.711	0.369	3.052
Trait Number_P	0.097	-0.484	0.677

3.5 – Meta-analysis of Angle Differences Within Environments

First, we can have a look at the dataset more carefully before running the analysis.

    kable(angle_data_within) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

X	study	env	ES	S_var	environment	authors	year	avg_n_env	trait_num	genus	species	design	group	env.comp	species_full	stdy	err	wi
1	RR073	B	-0.3419286	0.5873615	E1	Evans et al.	2015	114.5	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	1	1.7025290
2	RR073	S	-0.3158065	0.0463421	E2	Evans et al.	2015	105.0	8	Poecilia	reticulata	half-sib	fish	BS	Poecilia_reticulata	RR073	2	21.5786463
3	RR076	B	-2.9218833	0.0317398	E1	Guan et al.	2016	1097.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	3	31.5061726
4	RR076	S	-1.9911908	0.0618567	E2	Guan et al.	2016	1497.0	3	Scophthalmus	maximus	half-sib	fish	BS	Scophthalmus_maximus	RR076	4	16.1663968
5	RR087a	A	-0.9597197	0.1376774	E1	Sae-Lim et al.	2013	1232.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	5	7.2633537
6	RR087a	N	-1.4588856	0.0408052	E2	Sae-Lim et al.	2013	995.5	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	6	24.5066969
7	RR087b	A	-0.9584868	0.1239307	E1	Sae-Lim et al.	2013	1232.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	7	8.0690287
8	RR087b	N	-0.8841868	0.0427569	E2	Sae-Lim et al.	2013	1445.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	8	23.3880215
9	RR087c	A	-0.9148208	0.1567327	E1	Sae-Lim et al.	2013	1232.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	9	6.3802888
10	RR087c	N	-2.0659812	0.8070163	E2	Sae-Lim et al.	2013	959.0	3	Oncorhynchus	mykiss	full-sib	fish	AN	Oncorhynchus_mykiss	RR087	10	1.2391324
11	sR097	B	-1.6290017	0.1391402	E1	Kasule	1991	169.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	11	7.1869948
12	sR097	S	-2.8359926	1.3165859	E2	Kasule	1991	169.0	2	Dysdercus	fasciatus	half-sib	insect	BS	Dysdercus_fasciatus	sR097	12	0.7595403
13	sR110b	B	-3.2403778	0.7802790	E1	Blanckenhorn & Heyland	2004	82.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	13	1.2815928
14	sR110b	S	-2.9282182	1.4052676	E2	Blanckenhorn & Heyland	2004	62.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	14	0.7116082
15	sR110c	B	-3.7627363	1.2093555	E1	Blanckenhorn & Heyland	2004	105.5	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	15	0.8268867
16	sR110c	S	-3.2047958	1.3381308	E2	Blanckenhorn & Heyland	2004	81.0	2	Scathophaga	stercoraria	half-sib	insect	BS	Scathophaga_stercoraria	sR110	16	0.7473111
17	sR115	A	-2.4518216	0.2393191	E1	Conner et al.	2003	741.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	17	4.1785219
18	sR115	N	-0.1966421	0.0056296	E2	Conner et al.	2003	740.5	6	Raphanus	raphanistrum	half-sib	plant	AN	Raphanus_raphanistrum	sR115	18	177.6331442
19	sR116	A	-2.0465360	0.0292682	E1	Delcourt et al.	2010	1840.5	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	19	34.1667396
20	sR116	N	-1.8720304	0.0413366	E2	Delcourt et al.	2010	1824.0	8	Drosophila	serrata	half-sib	insect	AN	Drosophila_serrata	sR116	20	24.1916409
21	sR145	B	0.0908591	0.7744759	E1	Tucic et al.	1991	101.0	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	21	1.2911958
22	sR145	S	0.7859399	0.5287318	E2	Tucic et al.	1991	83.5	8	Acanthoscelides	obtectus	full-sib	insect	BS	Acanthoscelides_obtectus	sR145	22	1.8913179
23	WB002a	A	-1.8890462	0.2480546	E1	Begin & Roff	2001	372.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	23	4.0313701
24	WB002a	N	-1.8143146	0.1651411	E2	Begin & Roff	2001	272.0	5	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB002	24	6.0554263
25	WB003	A	-2.9213749	0.1765009	E1	Begin et al.	2004	294.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	25	5.6656937
26	WB003	N	-3.0281529	0.2354223	E2	Begin et al.	2004	387.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	26	4.2476864
27	WB003_2	A	-2.9298980	0.1733582	E1	Begin et al.	2004	294.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	27	5.7684033
28	WB003_2	N	-2.8386106	0.2638060	E2	Begin et al.	2004	462.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	28	3.7906645
29	WB003a	A	-2.6253466	0.1477390	E1	Begin et al.	2004	153.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	29	6.7686951
30	WB003a	N	-2.5919519	0.1156100	E2	Begin et al.	2004	82.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	30	8.6497698
31	WB003a_2	A	-2.6548252	0.1429089	E1	Begin et al.	2004	153.5	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	31	6.9974646
32	WB003a_2	N	-3.0880015	0.1029547	E2	Begin et al.	2004	178.0	5	Gryllus	firmus	full-sib	insect	AN	Gryllus_firmus	WB003	32	9.7130125
33	WB007	B	1.2753543	0.0083553	E1	Brock et al.	2010	366.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	33	119.6844085
34	WB007	S	2.4047661	0.1121431	E2	Brock et al.	2010	327.5	11	Brassica	rapa	half-sib	plant	BS	Brassica_rapa	WB007	34	8.9171751
35	WB010	A	5.5349758	1.2283841	E1	Collins et al.	1999	840.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	35	0.8140776
36	WB010	N	3.2749170	0.3578075	E2	Collins et al.	1999	252.5	4	Achroia	grisella	half-sib	insect	AN	Achroia_grisella	WB010	36	2.7947988
37	WB011	B	2.7850909	0.0260469	E1	Czesak & Fox	2003	1535.5	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	37	38.3922801
38	WB011	S	1.6471936	0.0029787	E2	Czesak & Fox	2003	1535.5	2	Stator	limbatus	half-sib	insect	BS	Stator_limbatus	WB011	38	335.7189771
39	WB028	B	-2.3077811	0.6592608	E1	Klause & Morin	2001	330.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	39	1.5168505
40	WB028	S	-4.5762373	1.0977522	E2	Klause & Morin	2001	330.0	2	Priophorus	pallipes	full-sib	insect	BS	Priophorus_pallipes	WB028	40	0.9109524
41	WB029	B	0.7016899	0.0028168	E1	King et al.	2011	862.5	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	41	355.0121381
42	WB029	S	-0.6077572	0.0018738	E2	King et al.	2011	895.5	3	Gryllus	firmus	half-sib	insect	BS	Gryllus_firmus	WB029	42	533.6691304
43	WB033	A	-4.6308526	0.0834496	E1	Lau et al.	2014	293.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	43	11.9832766
44	WB033	N	-5.0518975	0.2464121	E2	Lau et al.	2014	280.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	44	4.0582423
45	WB033a	A	-4.1964803	0.0220943	E1	Lau et al.	2014	259.5	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	45	45.2604759
46	WB033a	N	-5.2362594	0.1008700	E2	Lau et al.	2014	265.0	4	Arabidopsis	thaliana	half-sib	plant	AN	Arabidopsis_thaliana	WB033	46	9.9137479
47	WB035	B	-2.8244414	0.0259261	E1	Messina et al.	2013	1316.5	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	47	38.5711642
48	WB035	S	-1.0705922	0.0404241	E2	Messina et al.	2013	1319.5	3	Callosobruchus	maculatus	half-sib	insect	BS	Callosobruchus_maculatus	WB035	48	24.7377434
49	WB038	B	-3.0948590	0.2109398	E1	Rauter&Moore	2002	387.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	49	4.7406899
50	WB038	S	-2.4697815	0.1683493	E2	Rauter&Moore	2002	387.0	4	Nicrophorus	pustulatus	half-sib	insect	BS	Nicrophorus _pustulatus	WB038	50	5.9400312
51	WB045	A	-2.4725428	0.1708528	E1	Simons & Roff	1996	209.5	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	51	5.8529902
52	WB045	N	-0.2165210	0.0005315	E2	Simons & Roff	1996	200.5	7	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	52	1881.2980685
53	WB045a	A	-2.2064629	0.2477518	E1	Simons & Roff	1996	215.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	53	4.0362980
54	WB045a	N	-3.0483883	0.2677812	E2	Simons & Roff	1996	205.0	6	Gryllus	pennsylvanicus	full-sib	insect	AN	Gryllus_pennsylvanicus	WB045	54	3.7343921
55	WB053	A	1.1219771	0.0184589	E1	Via & Conner	1995	252.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	55	54.1744813
56	WB053	N	0.6552539	0.0227371	E2	Via & Conner	1995	252.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	56	43.9809635
57	WB053_3	A	0.8219987	0.0188514	E1	Via & Conner	1995	252.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	57	53.0465318
58	WB053_3	N	1.7625772	0.7594263	E2	Via & Conner	1995	252.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	58	1.3167835
59	WB053b	A	0.6628335	0.0030724	E1	Via & Conner	1995	280.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	59	325.4826934
60	WB053b	N	1.2106564	0.0342508	E2	Via & Conner	1995	280.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	60	29.1963872
61	WB053b_3	A	0.4629305	0.0030871	E1	Via & Conner	1995	280.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	61	323.9316315
62	WB053b_3	N	0.4594170	0.0387154	E2	Via & Conner	1995	280.0	2	Tribolium	castaneum	half-sib	insect	AN	Tribolium_castaneum	WB053	62	25.8295272
63	WB057	A	-2.7696908	0.1094535	E1	Grill et al.	1997	149.5	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	63	9.1362990
64	WB057	N	-2.4674966	0.3017626	E2	Grill et al.	1997	151.5	4	Harmonia	axyridis	full-sib	insect	AN	Harmonia_axyridis	WB057	64	3.3138633

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

     model_intcp_angle_within_het <- rma.mv(ES ~ 1, V = angle_data$S_var,   random = list(~1|stdy, ~1|err), data = angle_data_within)

            stdy_sig_angle_within <- model_intcp_angle_within_het$sigma2[1]
             err_sig_angle_within <- model_intcp_angle_within_het$sigma2[2]

        # Calculate total sampling variance
            wi_angle_within <- 1 / angle_data_within$S_var
            vi_angle_within <- sum(wi_angle_within) * ( nrow(angle_data_within) - 1) / ( (sum(wi_angle_within)^2) - sum(wi_angle_within))

        # Heterogeneity
    I2total_angle_within <- round((stdy_sig_angle_within + err_sig_angle_within) / (stdy_sig_angle_within + err_sig_angle_within + vi_angle_within), digits = 2)

    I2stdy_angle_within <- round((stdy_sig_angle_within) / (stdy_sig_angle_within + err_sig_angle_within + vi_angle_within), digits = 2)

Again, this will provide estimates of heterogeneity as the I2total_angle_within = 1 and between study heterogeneity as I2stdy_angle_within = 0.9. While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the environment type of matrix (non-novel or novel), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

    # Angle within environment multi-level meta regression models. Here we want to account for the fact that we have two environments per study and so these are correlated. 
        model_Angle_withinenv <- rma.mv(ES ~ environment + design + env.comp + scale(trait_num), V = angle_data_within$S_var, random = list(~1 + environment|stdy, ~1|err), data = angle_data_within)
            model_Angle_withinenv_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = angle_data_within$S_var, random = list(~1|stdy, ~1|err), data = angle_data_within)
          model_Angle_withinenv_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = angle_data_within$S_var, random = list(~1|stdy, ~1|err), data = angle_data_within)

        n_AN_S_angle_within <- table(angle_data_within$env.comp)
        n_design_angle_within <- table(angle_data_within$design)

        # Marginal means
        angle_within_G_AN <- marginal(model_Angle_withinenv, "design", data = angle_data_within)
        angle_within_G_S  <- marginal(model_Angle_withinenv_stress, "design", data = angle_data_within)

        env_cmp_angle_within <- rbind(angle_within_G_AN, angle_within_G_S)
        env_cmp_angle_within <- cbind(env_cmp_angle_within, n_AN_S_angle_within)

        Angle_within_G_full <- marginal(model_Angle_withinenv, "env.comp", angle_data_within)
        Angle_within_G_half <- marginal(model_Angle_withinenv_half, "env.comp", angle_data_within)
          
        design_angle_within <- rbind(Angle_within_G_full, Angle_within_G_half)
        design_angle_within <- cbind(design_angle_within, n_design_angle_within)

        angle_within_marg_means <- rbind(env_cmp_angle_within, design_angle_within)

Once we have these models run we can have a look at the table of estimates:

    table_angle_within <- data.frame(Est. = round(model_Angle_withinenv$b, digits =3), L95CI = round(model_Angle_withinenv$ci.lb, digits=3), U95CI = round(model_Angle_withinenv$ci.ub, digits =3))
    rownames(table_angle_within) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Stressful Environment", "Trait Number")

    kable(table_angle_within) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept	-2.130	-3.858	-0.402
Novel Environment	0.147	-0.305	0.599
Half-sib Breeding Design	1.078	-1.024	3.180
Stressful Environment	0.095	-1.886	2.076
Trait Number	0.754	-0.050	1.558

4.0 – Meta-analysis 2: Evolvability of environmentally induced phenotypes

4.1 – Load Data

              evolve_max_data <- read.csv("./data_evolvability/evolve_max_data_standardised.csv")
         evolve_max_data_test <- read.csv("./data_evolvability/evolve_max_data_test.csv")
            evolve_angle_data <- read.csv("./data_evolvability/evolve_angle_data_standardised.csv")

            evolve_max_data$stress <- ifelse(evolve_max_data$env == "S", "S", "A")
            evolve_angle_data$stress <- ifelse(evolve_angle_data$env == "S", "S", "A")

4.2 – Meta-analysis of \(\pi_{e}\)

Again, we can view the raw data as normal.

    kable(evolve_max_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

authors	year	study	env	type	n	trait_num	genus	species	design	group	name.type	name.env	ES	S_var	err	environment	stdy
Beacham	1988	RR062	N	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062.G	RR062.N	0.5428624	0.1240245	1	Novel	RR062
Beacham	1988	RR062_b	N	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062_b.G	RR062_b.N	0.2101515	1.7883816	2	Novel	RR062
Beacham	1988	RR062b	N	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b.G	RR062b.N	-1.2101635	0.8238148	3	Novel	RR062
Beacham	1988	RR062b_b	N	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b_b.G	RR062b_b.N	-4.5066487	0.2344635	4	Novel	RR062
Evans et al.	2015	RR073	S	G	64	8	Poecilia	reticulata	half-sib	fish	RR073.G	RR073.S	0.4505784	0.0172563	5	Novel	RR073
Guan et al.	2016	RR076	S	G	69	3	Scophthalmus	maximus	half-sib	fish	RR076.G	RR076.S	1.3571034	0.0152012	6	Novel	RR076
Sae-Lim et al.	2013	RR087a	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087a.G	RR087a.N	0.6551657	0.0294507	7	Novel	RR087
Sae-Lim et al.	2013	RR087b	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087b.G	RR087b.N	-0.1804580	0.0015965	8	Novel	RR087
Sae-Lim et al.	2013	RR087c	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087c.G	RR087c.N	0.6335908	0.1395597	9	Novel	RR087
Kasule	1991	sR097	S	G	26	2	Dysdercus	fasciatus	half-sib	insect	sR097.G	sR097.S	5.0286449	4.4480211	10	Novel	sR097
Blanckenhorn & Heyland	2004	sR110b	S	G	27	2	Scathophaga	stercoraria	half-sib	insect	sR110b.G	sR110b.S	1.7989274	0.6511523	11	Novel	sR110
Blanckenhorn & Heyland	2004	sR110c	S	G	33	2	Scathophaga	stercoraria	half-sib	insect	sR110c.G	sR110c.S	1.4345503	0.2583591	12	Novel	sR110
Conner et al.	2003	sR115	N	G	593	6	Raphanus	raphanistrum	half-sib	plant	sR115.G	sR115.N	0.8059430	0.0037772	13	Novel	sR115
Delcourt et al.	2010	sR116	N	G	365	8	Drosophila	serrata	half-sib	insect	sR116.G	sR116.N	-4.0826705	0.0077509	14	Novel	sR116
Tucic et al.	1991	sR145	S	G	43	8	Acanthoscelides	obtectus	full-sib	insect	sR145.G	sR145.S	0.0775910	0.0771394	15	Novel	sR145
Auld	2010	WB001	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001.G	WB001.S	0.2444737	0.6876215	16	Novel	WB001
Auld	2010	WB001_b	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001_b.G	WB001_b.S	2.5420243	0.7899251	17	Novel	WB001
Auld	2010	WB001_c	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001_c.G	WB001_c.S	1.3910750	0.1019580	18	Novel	WB001
Begin & Roff	2001	WB002a	N	G	39	5	Gryllus	pennsylvanicus	full-sib	insect	WB002a.G	WB002a.N	0.9910524	0.0537902	19	Novel	WB002
Begin et al.	2004	WB003	N	G	62	5	Gryllus	firmus	full-sib	insect	WB003.G	WB003.N	-3.0566626	0.0450678	20	Novel	WB003
Begin et al.	2004	WB003_2	N	G	60	5	Gryllus	firmus	full-sib	insect	WB003_2.G	WB003_2.N	-2.3639847	0.0373027	21	Novel	WB003
Begin et al.	2004	WB003a	N	G	41	5	Gryllus	firmus	full-sib	insect	WB003a.G	WB003a.N	0.8568436	0.0204842	22	Novel	WB003
Begin et al.	2004	WB003a_2	N	G	32	5	Gryllus	firmus	full-sib	insect	WB003a_2.G	WB003a_2.N	-0.4183020	0.0048221	23	Novel	WB003
Brock et al.	2010	WB007	S	G	131	11	Brassica	rapa	half-sib	plant	WB007.G	WB007.S	-3.1248238	0.0085661	24	Novel	WB007
Collins et al.	1999	WB010	N	G	62	4	Achroia	grisella	half-sib	insect	WB010.G	WB010.N	1.1622610	0.0075125	25	Novel	WB010
Czesak & Fox	2003	WB011	S	G	404	2	Stator	limbatus	half-sib	insect	WB011.G	WB011.S	-2.9885446	0.0001467	26	Novel	WB011
Delcourt et al.	2009	WB014	N	G	317	2	Drosophila	serrata	half-sib	insect	WB014.G	WB014.N	5.2775656	0.0003394	27	Novel	WB014
Engqvist	2007	WB016	S	G	27	2	Panorpa	cognata	full-sib	insect	WB016.G	WB016.S	10.0864227	2.2858387	28	Novel	WB016
Guntrip et al.	1997	WB022	N	G	40	2	Callosobruchus	maculatus	half-sib	insect	WB022.G	WB022.N	0.8558881	0.1851336	29	Novel	WB022
Holloway et al.	1990	WB025	N	G	69	4	Sitophilus	oryzae	half-sib	insect	WB025.G	WB025.N	1.2765941	0.0699061	30	Novel	WB025
Klause & Morin	2001	WB028	S	G	60	2	Priophorus	pallipes	full-sib	insect	WB028.G	WB028.S	-1.1091105	0.0060447	31	Novel	WB028
King et al.	2011	WB029	S	G	63	3	Gryllus	firmus	half-sib	insect	WB029.G	WB029.S	-2.9537509	0.0131215	32	Novel	WB029
Lau et al.	2014	WB033	N	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033.G	WB033.N	1.0437744	0.0000261	33	Novel	WB033
Lau et al.	2014	WB033a	N	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033a.G	WB033a.N	-2.9124026	0.0000010	34	Novel	WB033
Messina et al.	2013	WB035	S	G	432	3	Callosobruchus	maculatus	half-sib	insect	WB035.G	WB035.S	5.5948655	0.7656785	35	Novel	WB035
Paccard et al.	2013	WB036	N	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036.G	WB036.N	0.7126878	0.0132725	36	Novel	WB036
Paccard et al.	2013	WB036A	N	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036A.G	WB036A.N	1.8722282	0.0597320	37	Novel	WB036
Rauter&Moore	2002	WB038	S	G	135	4	Nicrophorus	pustulatus	half-sib	insect	WB038.G	WB038.S	3.4314013	0.0984880	38	Novel	WB038
Relyea	2005	WB039	S	G	21	9	Rana	sylvatica	half-sib	amphibian	WB039.G	WB039.S	-0.1582131	0.1651519	39	Novel	WB039
Sherrard et al.	2009	WB044	S	G	26	10	Avena	barbata	half-sib	plant	WB044.G	WB044.S	-3.8866656	0.1485572	40	Novel	WB044
Simons & Roff	1996	WB045	N	G	39	7	Gryllus	pennsylvanicus	full-sib	insect	WB045.G	WB045.N	0.5875816	0.0024073	41	Novel	WB045
Simons & Roff	1996	WB045a	N	G	39	6	Gryllus	pennsylvanicus	full-sib	insect	WB045a.G	WB045a.N	3.2819619	0.0732472	42	Novel	WB045
Simonsen & Stinchcombe	2001	WB046	S	G	50	3	Ipomoea	hederacea	half-sib	plant	WB046.G	WB046.S	4.1499697	0.9443279	43	Novel	WB046
Via & Conner	1995	WB053	N	G	24	2	Tribolium	castaneum	half-sib	insect	WB053.G	WB053.N	0.6486913	0.0301427	44	Novel	WB053
Via & Conner	1995	WB053_3	N	G	24	2	Tribolium	castaneum	half-sib	insect	WB053_3.G	WB053_3.N	0.0056850	0.2296783	45	Novel	WB053
Via & Conner	1995	WB053b	N	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b.G	WB053b.N	0.7647140	0.0100357	46	Novel	WB053
Via & Conner	1995	WB053b_3	N	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b_3.G	WB053b_3.N	-0.5880988	0.0682137	47	Novel	WB053
Service & Rose	1985	WB056	N	G	64	2	Drosophila	melanogaster	half-sib	insect	WB056.G	WB056.N	1.8715821	0.4847043	48	Novel	WB056
Grill et al.	1997	WB057	N	G	26	4	Harmonia	axyridis	full-sib	insect	WB057.G	WB057.N	-1.3737456	0.0503796	49	Novel	WB057
Beacham	1988	RR062	A	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062.G	RR062.A	-1.8588986	0.4660717	50	Non-novel	RR062
Beacham	1988	RR062_b	A	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062_b.G	RR062_b.A	-0.4636018	0.6964657	51	Non-novel	RR062
Beacham	1988	RR062b	A	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b.G	RR062b.A	-2.4161392	0.4171573	52	Non-novel	RR062
Beacham	1988	RR062b_b	A	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b_b.G	RR062b_b.A	-0.2500228	0.5239350	53	Non-novel	RR062
Evans et al.	2015	RR073	B	G	68	8	Poecilia	reticulata	half-sib	fish	RR073.G	RR073.B	1.1675329	0.0670842	54	Non-novel	RR073
Guan et al.	2016	RR076	B	G	69	3	Scophthalmus	maximus	half-sib	fish	RR076.G	RR076.B	1.0152374	0.0082485	55	Non-novel	RR076
Sae-Lim et al.	2013	RR087a	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087a.G	RR087a.A	1.3459910	0.1654308	56	Non-novel	RR087
Sae-Lim et al.	2013	RR087b	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087b.G	RR087b.A	-0.0889492	0.0226032	57	Non-novel	RR087
Sae-Lim et al.	2013	RR087c	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087c.G	RR087c.A	0.5866717	0.1139000	58	Non-novel	RR087
Kasule	1991	sR097	B	G	26	2	Dysdercus	fasciatus	half-sib	insect	sR097.G	sR097.B	4.3263515	1.7230775	59	Non-novel	sR097
Blanckenhorn & Heyland	2004	sR110b	B	G	37	2	Scathophaga	stercoraria	half-sib	insect	sR110b.G	sR110b.B	-0.1226381	0.0067473	60	Non-novel	sR110
Blanckenhorn & Heyland	2004	sR110c	B	G	40	2	Scathophaga	stercoraria	half-sib	insect	sR110c.G	sR110c.B	0.6388066	0.0367883	61	Non-novel	sR110
Conner et al.	2003	sR115	A	G	350	6	Raphanus	raphanistrum	half-sib	plant	sR115.G	sR115.A	3.4973585	0.1625815	62	Non-novel	sR115
Delcourt et al.	2010	sR116	A	G	366	8	Drosophila	serrata	half-sib	insect	sR116.G	sR116.A	-4.4340787	0.0079699	63	Non-novel	sR116
Tucic et al.	1991	sR145	B	G	52	8	Acanthoscelides	obtectus	full-sib	insect	sR145.G	sR145.B	0.3300179	0.1013896	64	Non-novel	sR145
Auld	2010	WB001	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001.G	WB001.B	-1.2919352	0.2238658	65	Non-novel	WB001
Auld	2010	WB001_b	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001_b.G	WB001_b.B	1.9766137	0.7477113	66	Non-novel	WB001
Auld	2010	WB001_c	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001_c.G	WB001_c.B	1.1365886	0.6238066	67	Non-novel	WB001
Begin & Roff	2001	WB002a	A	G	39	5	Gryllus	pennsylvanicus	full-sib	insect	WB002a.G	WB002a.A	0.9551444	0.0507751	68	Non-novel	WB002
Begin et al.	2004	WB003	A	G	56	5	Gryllus	firmus	full-sib	insect	WB003.G	WB003.A	-3.0254976	0.0642904	69	Non-novel	WB003
Begin et al.	2004	WB003_2	A	G	56	5	Gryllus	firmus	full-sib	insect	WB003_2.G	WB003_2.A	-2.5292244	0.0540068	70	Non-novel	WB003
Begin et al.	2004	WB003a	A	G	45	5	Gryllus	firmus	full-sib	insect	WB003a.G	WB003a.A	0.4031230	0.0321589	71	Non-novel	WB003
Begin et al.	2004	WB003a_2	A	G	45	5	Gryllus	firmus	full-sib	insect	WB003a_2.G	WB003a_2.A	-0.5568401	0.0224074	72	Non-novel	WB003
Brock et al.	2010	WB007	B	G	132	11	Brassica	rapa	half-sib	plant	WB007.G	WB007.B	-0.9073448	0.0187114	73	Non-novel	WB007
Collins et al.	1999	WB010	A	G	141	4	Achroia	grisella	half-sib	insect	WB010.G	WB010.A	1.4115091	0.0019616	74	Non-novel	WB010
Czesak & Fox	2003	WB011	B	G	404	2	Stator	limbatus	half-sib	insect	WB011.G	WB011.B	-3.0142257	0.0006719	75	Non-novel	WB011
Delcourt et al.	2009	WB014	A	G	319	2	Drosophila	serrata	half-sib	insect	WB014.G	WB014.A	5.2246319	0.0017553	76	Non-novel	WB014
Engqvist	2007	WB016	B	G	27	2	Panorpa	cognata	full-sib	insect	WB016.G	WB016.B	12.5280490	4.9242628	77	Non-novel	WB016
Guntrip et al.	1997	WB022	A	G	40	2	Callosobruchus	maculatus	half-sib	insect	WB022.G	WB022.A	-1.3775179	0.0789878	78	Non-novel	WB022
Holloway et al.	1990	WB025	A	G	65	4	Sitophilus	oryzae	half-sib	insect	WB025.G	WB025.A	2.1329860	0.0104086	79	Non-novel	WB025
Klause & Morin	2001	WB028	B	G	60	2	Priophorus	pallipes	full-sib	insect	WB028.G	WB028.B	-1.2711403	0.0982403	80	Non-novel	WB028
King et al.	2011	WB029	B	G	63	3	Gryllus	firmus	half-sib	insect	WB029.G	WB029.B	-2.7665833	0.0132234	81	Non-novel	WB029
Lau et al.	2014	WB033	A	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033.G	WB033.A	1.0449722	0.0000227	82	Non-novel	WB033
Lau et al.	2014	WB033a	A	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033a.G	WB033a.A	-2.9069230	0.0000013	83	Non-novel	WB033
Messina et al.	2013	WB035	B	G	432	3	Callosobruchus	maculatus	half-sib	insect	WB035.G	WB035.B	0.9121341	0.0020951	84	Non-novel	WB035
Paccard et al.	2013	WB036	A	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036.G	WB036.A	0.7286765	0.1116486	85	Non-novel	WB036
Paccard et al.	2013	WB036A	A	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036A.G	WB036A.A	2.3949669	0.0263805	86	Non-novel	WB036
Rauter&Moore	2002	WB038	B	G	135	4	Nicrophorus	pustulatus	half-sib	insect	WB038.G	WB038.B	3.7426764	0.1518899	87	Non-novel	WB038
Relyea	2005	WB039	B	G	21	9	Rana	sylvatica	half-sib	amphibian	WB039.G	WB039.B	-0.3485422	0.1641508	88	Non-novel	WB039
Sherrard et al.	2009	WB044	B	G	26	10	Avena	barbata	half-sib	plant	WB044.G	WB044.B	-4.1224059	0.1631106	89	Non-novel	WB044
Simons & Roff	1996	WB045	A	G	39	7	Gryllus	pennsylvanicus	full-sib	insect	WB045.G	WB045.A	0.1642072	0.0270893	90	Non-novel	WB045
Simons & Roff	1996	WB045a	A	G	39	6	Gryllus	pennsylvanicus	full-sib	insect	WB045a.G	WB045a.A	3.2745521	0.2456119	91	Non-novel	WB045
Simonsen & Stinchcombe	2001	WB046	B	G	50	3	Ipomoea	hederacea	half-sib	plant	WB046.G	WB046.B	2.0149730	0.2862508	92	Non-novel	WB046
Via & Conner	1995	WB053	A	G	24	2	Tribolium	castaneum	half-sib	insect	WB053.G	WB053.A	1.0237109	0.0023335	93	Non-novel	WB053
Via & Conner	1995	WB053_3	A	G	24	2	Tribolium	castaneum	half-sib	insect	WB053_3.G	WB053_3.A	-1.1432159	0.0018889	94	Non-novel	WB053
Via & Conner	1995	WB053b	A	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b.G	WB053b.A	1.0418965	0.0006471	95	Non-novel	WB053
Via & Conner	1995	WB053b_3	A	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b_3.G	WB053b_3.A	-1.1265929	0.0005348	96	Non-novel	WB053
Service & Rose	1985	WB056	A	G	64	2	Drosophila	melanogaster	half-sib	insect	WB056.G	WB056.A	3.4838459	0.0562978	97	Non-novel	WB056
Grill et al.	1997	WB057	A	G	22	4	Harmonia	axyridis	full-sib	insect	WB057.G	WB057.A	-0.6138681	0.0260655	98	Non-novel	WB057

We can now estimate heterogeneity using a multi-level meta-analytic model

    mod_evolve_max_int <-  rma.mv(ES ~ 1, V = evolve_max_data$S_var, random = list(~1|stdy, ~1|err), data = evolve_max_data)

            stdy_evolve_max_int <- mod_evolve_max_int$sigma2[1]
             err_evolve_max_int <- mod_evolve_max_int$sigma2[2]
    
    # Calculate total sampling variance
            wi_evolve_max <- 1 / evolve_max_data$S_var
            vi_evolve_max <- sum(wi_evolve_max) * ( nrow(evolve_max_data) - 1) / ( (sum(wi_evolve_max)^2) - sum(wi_evolve_max))

            I2total_evolve_max_data <- round((stdy_evolve_max_int + err_evolve_max_int) / (stdy_evolve_max_int + err_evolve_max_int + vi_evolve_max), digits = 2)

            I2stdy_evolve_max_data <- round((stdy_evolve_max_int) / (stdy_evolve_max_int + err_evolve_max_int + vi_evolve_max), digits = 2)

OK, but again, we have a priori moderators we think should explain these

            mod_evolve_max_trait_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)
            mod_evolve_max_trait_env_stress <-  rma.mv(ES ~ 1 + environment + design + stress + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)

        # calculate marginal mean for non-novel, marginalised over design type. First, refit model
            mod_evolve_max_trait_env_novel <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)          
            
            marinal_means_evolve_max_A_N <- rbind(marginal(mod_evolve_max_trait_env, "design", data = evolve_max_data), marginal(mod_evolve_max_trait_env_novel, "design", data = evolve_max_data))
            n_A_N <- table(evolve_max_data$environment)
            marinal_means_evolve_max_A_N <- cbind(marinal_means_evolve_max_A_N, n_A_N)
            
        # Now we can do this over design
            mod_evolve_max_trait_env_half <-  rma.mv(ES ~ 1 + environment + relevel(design, ref = "half-sib") + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)
            marinal_means_evolve_max_design <- rbind(marginal(mod_evolve_max_trait_env, "environment", data = evolve_max_data), marginal(mod_evolve_max_trait_env_half, "environment", data = evolve_max_data))
            n_design <- table(evolve_max_data$design)
            marinal_means_evolve_max_design <- cbind(marinal_means_evolve_max_design, n_design)

            marg_evolve_max <- rbind(marinal_means_evolve_max_A_N, marinal_means_evolve_max_design)

Once we have these models run we can have a look at the table of estimates:

    table_max_evolve <- data.frame(Est. = round(mod_evolve_max_trait_env$b, digits =3), L95CI = round(mod_evolve_max_trait_env$ci.lb, digits=3), U95CI = round(mod_evolve_max_trait_env$ci.ub, digits =3))
    rownames(table_max_evolve) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")   
    kable(table_max_evolve) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept	1.268	-0.356	2.893
Novel Environment	0.092	-0.417	0.601
Half-sib Breeding Design	-0.754	-2.685	1.176
Trait Number	-1.073	-1.873	-0.274

To test whether there is more additive genetic variation along the plasticity vector than expected we also run the model for the difference between \(\pi_{e}\) and \(\pi_{0}\).

        mod_evolve_max_test_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data_test$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data_test)

Once we have this model run we can have a look at the table of estimates:

4.3 – Table S1

Table S1 – Estimates testing whether the angle of gmax is aligned more than simply by chance. We are specifically interested in the intercept here and when the 95% confidence intervals exclude zero, there is evidence that alignment between plasticity vector and \(g_{max}\) is more then what we expect by chance alone.

    table_max_evolve_test <- data.frame(Est. = round(mod_evolve_max_test_env$b, digits =3), L95CI = round(mod_evolve_max_test_env$ci.lb, digits=3), U95CI = round(mod_evolve_max_test_env$ci.ub, digits =3))
    rownames(table_max_evolve_test) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")  
    kable(table_max_evolve_test) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept	1.108	0.393	1.822
Novel Environment	0.040	-0.241	0.321
Half-sib Breeding Design	-0.573	-1.420	0.274
Trait Number	0.164	-0.192	0.521

4.4 – Meta-analysis of \(\theta_{e}\)

In case a reader is interested, we’ll still spit out the raw data:

    kable(evolve_angle_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")

authors	year	study	env	type	n	trait_num	genus	species	design	group	name.type	name.env	ES	S_var	err	environment	stdy
Beacham	1988	RR062	N	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062.G	RR062.N	-0.2804839	0.0692202	1	Novel	RR062
Beacham	1988	RR062_b	N	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062_b.G	RR062_b.N	1.0696574	2.1866890	2	Novel	RR062
Beacham	1988	RR062b	N	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b.G	RR062b.N	2.1656787	1.4562036	3	Novel	RR062
Beacham	1988	RR062b_b	N	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b_b.G	RR062b_b.N	3.0168584	0.2128005	4	Novel	RR062
Evans et al.	2015	RR073	S	G	64	8	Poecilia	reticulata	half-sib	fish	RR073.G	RR073.S	-0.2136747	0.0122647	5	Novel	RR073
Guan et al.	2016	RR076	S	G	69	3	Scophthalmus	maximus	half-sib	fish	RR076.G	RR076.S	-0.8207908	0.0063809	6	Novel	RR076
Sae-Lim et al.	2013	RR087a	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087a.G	RR087a.N	-0.2667054	0.0127446	7	Novel	RR087
Sae-Lim et al.	2013	RR087b	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087b.G	RR087b.N	0.1299726	0.0008169	8	Novel	RR087
Sae-Lim et al.	2013	RR087c	N	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087c.G	RR087c.N	0.2813025	0.1212850	9	Novel	RR087
Kasule	1991	sR097	S	G	26	2	Dysdercus	fasciatus	half-sib	insect	sR097.G	sR097.S	-2.7552460	1.2302908	10	Novel	sR097
Blanckenhorn & Heyland	2004	sR110b	S	G	27	2	Scathophaga	stercoraria	half-sib	insect	sR110b.G	sR110b.S	-0.8894004	0.2351409	11	Novel	sR110
Blanckenhorn & Heyland	2004	sR110c	S	G	33	2	Scathophaga	stercoraria	half-sib	insect	sR110c.G	sR110c.S	-0.7073991	0.1014419	12	Novel	sR110
Conner et al.	2003	sR115	N	G	593	6	Raphanus	raphanistrum	half-sib	plant	sR115.G	sR115.N	-0.4036144	0.0043551	13	Novel	sR115
Delcourt et al.	2010	sR116	N	G	365	8	Drosophila	serrata	half-sib	insect	sR116.G	sR116.N	3.0267427	0.0073985	14	Novel	sR116
Tucic et al.	1991	sR145	S	G	43	8	Acanthoscelides	obtectus	full-sib	insect	sR145.G	sR145.S	0.2679014	0.1688566	15	Novel	sR145
Auld	2010	WB001	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001.G	WB001.S	1.6941891	1.8150091	16	Novel	WB001
Auld	2010	WB001_b	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001_b.G	WB001_b.S	-1.4002198	0.2552470	17	Novel	WB001
Auld	2010	WB001_c	S	G	30	3	Physa	acuta	full-sib	mollusc	WB001_c.G	WB001_c.S	-0.7901196	0.0400435	18	Novel	WB001
Begin & Roff	2001	WB002a	N	G	39	5	Gryllus	pennsylvanicus	full-sib	insect	WB002a.G	WB002a.N	-0.5388656	0.0264206	19	Novel	WB002
Begin et al.	2004	WB003	N	G	62	5	Gryllus	firmus	full-sib	insect	WB003.G	WB003.N	2.4367906	0.0410906	20	Novel	WB003
Begin et al.	2004	WB003_2	N	G	60	5	Gryllus	firmus	full-sib	insect	WB003_2.G	WB003_2.N	1.8226200	0.0232361	21	Novel	WB003
Begin et al.	2004	WB003a	N	G	41	5	Gryllus	firmus	full-sib	insect	WB003a.G	WB003a.N	-0.5049055	0.0085218	22	Novel	WB003
Begin et al.	2004	WB003a_2	N	G	32	5	Gryllus	firmus	full-sib	insect	WB003a_2.G	WB003a_2.N	0.2796241	0.0021071	23	Novel	WB003
Brock et al.	2010	WB007	S	G	131	11	Brassica	rapa	half-sib	plant	WB007.G	WB007.S	2.0510321	0.0038105	24	Novel	WB007
Collins et al.	1999	WB010	N	G	62	4	Achroia	grisella	half-sib	insect	WB010.G	WB010.N	-0.7145406	0.0036595	25	Novel	WB010
Czesak & Fox	2003	WB011	S	G	404	2	Stator	limbatus	half-sib	insect	WB011.G	WB011.S	1.8190447	0.0000475	26	Novel	WB011
Delcourt et al.	2009	WB014	N	G	317	2	Drosophila	serrata	half-sib	insect	WB014.G	WB014.N	-3.0455236	0.0000925	27	Novel	WB014
Engqvist	2007	WB016	S	G	27	2	Panorpa	cognata	full-sib	insect	WB016.G	WB016.S	-5.4896469	0.5746361	28	Novel	WB016
Guntrip et al.	1997	WB022	N	G	40	2	Callosobruchus	maculatus	half-sib	insect	WB022.G	WB022.N	-0.2539235	0.0864514	29	Novel	WB022
Holloway et al.	1990	WB025	N	G	69	4	Sitophilus	oryzae	half-sib	insect	WB025.G	WB025.N	-0.6432784	0.0496260	30	Novel	WB025
Klause & Morin	2001	WB028	S	G	60	2	Priophorus	pallipes	full-sib	insect	WB028.G	WB028.S	0.7417540	0.0024437	31	Novel	WB028
King et al.	2011	WB029	S	G	63	3	Gryllus	firmus	half-sib	insect	WB029.G	WB029.S	1.9110061	0.0050600	32	Novel	WB029
Lau et al.	2014	WB033	N	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033.G	WB033.N	-0.6590318	0.0000101	33	Novel	WB033
Lau et al.	2014	WB033a	N	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033a.G	WB033a.N	1.7678985	0.0000003	34	Novel	WB033
Messina et al.	2013	WB035	S	G	432	3	Callosobruchus	maculatus	half-sib	insect	WB035.G	WB035.S	-3.0488950	0.2429454	35	Novel	WB035
Paccard et al.	2013	WB036	N	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036.G	WB036.N	-0.4328543	0.0080883	36	Novel	WB036
Paccard et al.	2013	WB036A	N	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036A.G	WB036A.N	-1.1388280	0.0235012	37	Novel	WB036
Rauter&Moore	2002	WB038	S	G	135	4	Nicrophorus	pustulatus	half-sib	insect	WB038.G	WB038.S	-2.0022426	0.0314474	38	Novel	WB038
Relyea	2005	WB039	S	G	21	9	Rana	sylvatica	half-sib	amphibian	WB039.G	WB039.S	0.5733024	0.5624914	39	Novel	WB039
Sherrard et al.	2009	WB044	S	G	26	10	Avena	barbata	half-sib	plant	WB044.G	WB044.S	3.4875131	0.5825224	40	Novel	WB044
Simons & Roff	1996	WB045	N	G	39	7	Gryllus	pennsylvanicus	full-sib	insect	WB045.G	WB045.N	-0.3680302	0.0009672	41	Novel	WB045
Simons & Roff	1996	WB045a	N	G	39	6	Gryllus	pennsylvanicus	full-sib	insect	WB045a.G	WB045a.N	-1.9588013	0.0230411	42	Novel	WB045
Simonsen & Stinchcombe	2001	WB046	S	G	50	3	Ipomoea	hederacea	half-sib	plant	WB046.G	WB046.S	-2.3591948	0.2837434	43	Novel	WB046
Via & Conner	1995	WB053	N	G	24	2	Tribolium	castaneum	half-sib	insect	WB053.G	WB053.N	-0.3750967	0.0122579	44	Novel	WB053
Via & Conner	1995	WB053_3	N	G	24	2	Tribolium	castaneum	half-sib	insect	WB053_3.G	WB053_3.N	0.3731902	0.1313836	45	Novel	WB053
Via & Conner	1995	WB053b	N	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b.G	WB053b.N	-0.4465280	0.0040389	46	Novel	WB053
Via & Conner	1995	WB053b_3	N	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b_3.G	WB053b_3.N	0.9812257	0.0526150	47	Novel	WB053
Service & Rose	1985	WB056	N	G	64	2	Drosophila	melanogaster	half-sib	insect	WB056.G	WB056.N	-0.7985585	0.1836581	48	Novel	WB056
Grill et al.	1997	WB057	N	G	26	4	Harmonia	axyridis	full-sib	insect	WB057.G	WB057.N	1.0197897	0.0237647	49	Novel	WB057
Beacham	1988	RR062	A	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062.G	RR062.A	2.2968023	1.3286533	50	Non-novel	RR062
Beacham	1988	RR062_b	A	G	10	3	Oncorhynchus	gorbuscha	half-sib	fish	RR062_b.G	RR062_b.A	1.2979916	1.4852452	51	Non-novel	RR062
Beacham	1988	RR062b	A	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b.G	RR062b.A	2.4060164	0.8444542	52	Non-novel	RR062
Beacham	1988	RR062b_b	A	G	10	3	Oncorhynchus	keta	half-sib	fish	RR062b_b.G	RR062b_b.A	0.5874297	0.4801196	53	Non-novel	RR062
Evans et al.	2015	RR073	B	G	68	8	Poecilia	reticulata	half-sib	fish	RR073.G	RR073.B	-0.4779735	0.2397198	54	Non-novel	RR073
Guan et al.	2016	RR076	B	G	69	3	Scophthalmus	maximus	half-sib	fish	RR076.G	RR076.B	-0.6184308	0.0032563	55	Non-novel	RR076
Sae-Lim et al.	2013	RR087a	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087a.G	RR087a.A	-0.4300231	0.0797760	56	Non-novel	RR087
Sae-Lim et al.	2013	RR087b	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087b.G	RR087b.A	0.2516590	0.0208090	57	Non-novel	RR087
Sae-Lim et al.	2013	RR087c	A	G	100	3	Oncorhynchus	mykiss	full-sib	fish	RR087c.G	RR087c.A	0.1919137	0.0845678	58	Non-novel	RR087
Kasule	1991	sR097	B	G	26	2	Dysdercus	fasciatus	half-sib	insect	sR097.G	sR097.B	-2.4943512	0.4698627	59	Non-novel	sR097
Blanckenhorn & Heyland	2004	sR110b	B	G	37	2	Scathophaga	stercoraria	half-sib	insect	sR110b.G	sR110b.B	0.0981276	0.0027588	60	Non-novel	sR110
Blanckenhorn & Heyland	2004	sR110c	B	G	40	2	Scathophaga	stercoraria	half-sib	insect	sR110c.G	sR110c.B	-0.3332372	0.0152542	61	Non-novel	sR110
Conner et al.	2003	sR115	A	G	350	6	Raphanus	raphanistrum	half-sib	plant	sR115.G	sR115.A	-1.8909220	0.0593847	62	Non-novel	sR115
Delcourt et al.	2010	sR116	A	G	366	8	Drosophila	serrata	half-sib	insect	sR116.G	sR116.A	3.2396996	0.0071285	63	Non-novel	sR116
Tucic et al.	1991	sR145	B	G	52	8	Acanthoscelides	obtectus	full-sib	insect	sR145.G	sR145.B	0.4527517	0.8836976	64	Non-novel	sR145
Auld	2010	WB001	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001.G	WB001.B	2.6337736	1.1432721	65	Non-novel	WB001
Auld	2010	WB001_b	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001_b.G	WB001_b.B	-0.9856240	0.2666681	66	Non-novel	WB001
Auld	2010	WB001_c	B	G	30	3	Physa	acuta	full-sib	mollusc	WB001_c.G	WB001_c.B	-0.2218401	0.3625394	67	Non-novel	WB001
Begin & Roff	2001	WB002a	A	G	39	5	Gryllus	pennsylvanicus	full-sib	insect	WB002a.G	WB002a.A	-0.5221598	0.0250811	68	Non-novel	WB002
Begin et al.	2004	WB003	A	G	56	5	Gryllus	firmus	full-sib	insect	WB003.G	WB003.A	2.7580171	0.1227235	69	Non-novel	WB003
Begin et al.	2004	WB003_2	A	G	56	5	Gryllus	firmus	full-sib	insect	WB003_2.G	WB003_2.A	2.1443823	0.0515381	70	Non-novel	WB003
Begin et al.	2004	WB003a	A	G	45	5	Gryllus	firmus	full-sib	insect	WB003a.G	WB003a.A	-0.1690384	0.0148521	71	Non-novel	WB003
Begin et al.	2004	WB003a_2	A	G	45	5	Gryllus	firmus	full-sib	insect	WB003a_2.G	WB003a_2.A	0.4490378	0.0114060	72	Non-novel	WB003
Brock et al.	2010	WB007	B	G	132	11	Brassica	rapa	half-sib	plant	WB007.G	WB007.B	0.9015107	0.0265064	73	Non-novel	WB007
Collins et al.	1999	WB010	A	G	141	4	Achroia	grisella	half-sib	insect	WB010.G	WB010.A	-0.8776703	0.0007712	74	Non-novel	WB010
Czesak & Fox	2003	WB011	B	G	404	2	Stator	limbatus	half-sib	insect	WB011.G	WB011.B	1.8636743	0.0002255	75	Non-novel	WB011
Delcourt et al.	2009	WB014	A	G	319	2	Drosophila	serrata	half-sib	insect	WB014.G	WB014.A	-3.0175489	0.0004793	76	Non-novel	WB014
Engqvist	2007	WB016	B	G	27	2	Panorpa	cognata	full-sib	insect	WB016.G	WB016.B	-6.7136399	1.2335992	77	Non-novel	WB016
Guntrip et al.	1997	WB022	A	G	40	2	Callosobruchus	maculatus	half-sib	insect	WB022.G	WB022.A	1.3420900	0.0578258	78	Non-novel	WB022
Holloway et al.	1990	WB025	A	G	65	4	Sitophilus	oryzae	half-sib	insect	WB025.G	WB025.A	-1.3078633	0.0037445	79	Non-novel	WB025
Klause & Morin	2001	WB028	B	G	60	2	Priophorus	pallipes	full-sib	insect	WB028.G	WB028.B	2.1761862	0.4162618	80	Non-novel	WB028
King et al.	2011	WB029	B	G	63	3	Gryllus	firmus	half-sib	insect	WB029.G	WB029.B	1.8019485	0.0051812	81	Non-novel	WB029
Lau et al.	2014	WB033	A	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033.G	WB033.A	-0.6597876	0.0000088	82	Non-novel	WB033
Lau et al.	2014	WB033a	A	G	60	4	Arabidopsis	thaliana	half-sib	plant	WB033a.G	WB033a.A	1.7648037	0.0000004	83	Non-novel	WB033
Messina et al.	2013	WB035	B	G	432	3	Callosobruchus	maculatus	half-sib	insect	WB035.G	WB035.B	-0.5361305	0.0008421	84	Non-novel	WB035
Paccard et al.	2013	WB036	A	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036.G	WB036.A	-0.3174122	0.1071421	85	Non-novel	WB036
Paccard et al.	2013	WB036A	A	G	22	7	Arabidopsis	lyrata	full-sib	plant	WB036A.G	WB036A.A	-1.4644314	0.0091210	86	Non-novel	WB036
Rauter&Moore	2002	WB038	B	G	135	4	Nicrophorus	pustulatus	half-sib	insect	WB038.G	WB038.B	-2.1617375	0.0489554	87	Non-novel	WB038
Relyea	2005	WB039	B	G	21	9	Rana	sylvatica	half-sib	amphibian	WB039.G	WB039.B	0.7823934	0.6700404	88	Non-novel	WB039
Sherrard et al.	2009	WB044	B	G	26	10	Avena	barbata	half-sib	plant	WB044.G	WB044.B	4.1801808	1.0238365	89	Non-novel	WB044
Simons & Roff	1996	WB045	A	G	39	7	Gryllus	pennsylvanicus	full-sib	insect	WB045.G	WB045.A	-0.0302778	0.0128676	90	Non-novel	WB045
Simons & Roff	1996	WB045a	A	G	39	6	Gryllus	pennsylvanicus	full-sib	insect	WB045a.G	WB045a.A	-1.9201388	0.0779041	91	Non-novel	WB045
Simonsen & Stinchcombe	2001	WB046	B	G	50	3	Ipomoea	hederacea	half-sib	plant	WB046.G	WB046.B	-1.0651373	0.1040128	92	Non-novel	WB046
Via & Conner	1995	WB053	A	G	24	2	Tribolium	castaneum	half-sib	insect	WB053.G	WB053.A	-0.6445169	0.0009036	93	Non-novel	WB053
Via & Conner	1995	WB053_3	A	G	24	2	Tribolium	castaneum	half-sib	insect	WB053_3.G	WB053_3.A	0.7256621	0.0007313	94	Non-novel	WB053
Via & Conner	1995	WB053b	A	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b.G	WB053b.A	-0.6559237	0.0002504	95	Non-novel	WB053
Via & Conner	1995	WB053b_3	A	G	80	2	Tribolium	castaneum	half-sib	insect	WB053b_3.G	WB053b_3.A	0.7152326	0.0002055	96	Non-novel	WB053
Service & Rose	1985	WB056	A	G	64	2	Drosophila	melanogaster	half-sib	insect	WB056.G	WB056.A	-2.0719360	0.0170922	97	Non-novel	WB056
Grill et al.	1997	WB057	A	G	22	4	Harmonia	axyridis	full-sib	insect	WB057.G	WB057.A	0.4428487	0.0111591	98	Non-novel	WB057

We can now estimate heterogeneity using a multi-level meta-analytic model:

    mod_evolve_angle_int <-  rma.mv(ES ~ 1, V = evolve_angle_data$S_var, random = list(~1|stdy, ~1|err), data = evolve_angle_data)

        stdy_evolve_angle_int <- mod_evolve_angle_int$sigma2[1]
         err_evolve_angle_int <- mod_evolve_angle_int$sigma2[2]

            # Calculate total sampling variance
            wi_evolve_angle <- 1 / evolve_angle_data$S_var
            vi_volve_angle <- sum(wi_evolve_angle) * ( nrow(evolve_angle_data) - 1) / ( (sum(wi_evolve_angle)^2) - sum(wi_evolve_angle))

            I2total_evolve_max_data <- round((stdy_evolve_angle_int + err_sig_angle_within) / (stdy_evolve_angle_int + err_evolve_angle_int + vi_volve_angle), digits = 2)

            I2stdy_evolve_max_data <- round((stdy_evolve_angle_int) / (stdy_evolve_angle_int + err_evolve_angle_int + vi_volve_angle), digits = 2)

But again, we need to control for a number of factors likely to drive variation in effects.

    # Multi-level meta-regression models
        mod_evolve_angle_trait_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 + environment|stdy, ~1|err), data = evolve_angle_data)
        mod_evolve_angle_trait_env_stress <-  rma.mv(ES ~ 1 + environment + design + stress + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 + environment|stdy, ~1|err), data = evolve_angle_data)
        mod_evolve_angle_trait_env_stress_justNovel <-  rma.mv(ES ~ 1 + design + stress + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 |stdy, ~1|err), data = subset(evolve_angle_data, environment == "Novel"))

    # calculate marginal mean for non-novel, marginalised over design type. First, refit model
            mod_evolve_angle_trait_env_novel <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_angle_data)            
            
            marinal_means_evolve_angle_A_N <- rbind(marginal(mod_evolve_angle_trait_env, "design", data = evolve_angle_data), marginal(mod_evolve_angle_trait_env_novel, "design", data = evolve_angle_data))
            n_A_N <- table(evolve_angle_data$environment)
            marinal_means_evolve_angle_A_N <- cbind(marinal_means_evolve_angle_A_N, n_A_N)
            
        # Now we can do this over design
            mod_evolve_angle_trait_env_half <-  rma.mv(ES ~ 1 + environment + relevel(design, ref = "half-sib") + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_angle_data)
            marinal_means_evolve_angle_design <- rbind(marginal(mod_evolve_angle_trait_env, "environment", data = evolve_angle_data), marginal(mod_evolve_angle_trait_env_half, "environment", data = evolve_angle_data))
            n_design <- table(evolve_angle_data$design)
            marinal_means_evolve_angle_design <- cbind(marinal_means_evolve_angle_design, n_design)

            marg_evolve_angle <- rbind(marinal_means_evolve_angle_A_N, marinal_means_evolve_angle_design)

Once we have these models run we can have a look at the table of estimates:

    table_angle_evolve <- data.frame(Est. = round(mod_evolve_angle_trait_env$b, digits =3), L95CI = round(mod_evolve_angle_trait_env$ci.lb, digits=3), U95CI = round(mod_evolve_angle_trait_env$ci.ub, digits =3))
    rownames(table_angle_evolve) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")
    kable(table_angle_evolve) %>% kable_styling(font = 12, full_width = TRUE)

	Est.	L95CI	U95CI
Intercept	-0.522	-1.566	0.522
Novel Environment	-0.105	-0.461	0.251
Half-sib Breeding Design	0.450	-0.788	1.688
Trait Number	0.735	0.215	1.254

5.0 – Supplemental Figures

5.1 – Figure S1 – PRISMA

Figure S1 – PRISMA diagram for the systematic search.

5.2 – Figure S2 – Funnel plots

We explored publication bias using funnel plots, but there was no clear evidence across the effect sizes for publication bias:

Figure S2 – Funnel plots of sampling standard error and effect size. (A-C) ‘Gray’ points represent effect sizes for P, whereas ‘black’ points represent effect sizes for G. (D–F) ‘Orange’ points are effect sizes from novel environments, whereas ‘green’ points are effect sizes from non-novel environments. Note that for C) there was a significant change in angle across environments that was impacted also by the type of environment (stress or not). See main manuscript for more details.

5.3 – Code for Figure 2 - 5 – Main Paper

Figure 2

    #pdf(file = "Fig2.pdf", height = 7.055555, width = 17.513889)
    par(mfrow = c(1,2))

    # Just in case these slots are already in, lets just remove so there are no errors
    varalong_marg_means$name = NULL
    varalong_marg_means$obs = NULL

    marg_list_M1 <- list(SMD_marg_table, varalong_marg_means, angle_marg_means)
    marg_list_M1 <- lapply(marg_list_M1, function(x) cbind(x, name = c("Non-stressful", "Stressful", "Full-sib", "Half-sib")))
    marg_list_M1 <- lapply(marg_list_M1, function(x) cbind(x, obs = c(5,4,1,2)))

    labels = c( "Breeding Design", "Novel Environment")
    for(i in 1){
        par(mar = c(5,13,2,0))
        plot_marginal(marg_list_M1[[i]], "SDV", ylim = c(0,7), xexpan = c(-1.5, 1.5), cex.lab = 2, cex.axis = 1.5)
        #text(marg_list_M1[[i]]$Freq, x = 2, y = marg_list_M1[[i]]$obs, cex = 2)
        mtext(at = 3, side = 2, labels[1], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 6, side = 2, labels[2], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 7, side = 2, "A)", las = 1, font = 2.5, cex = 2.5, adj = 0.9, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Decreased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  "Total Variance", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1, expression(bold("Increased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, "Total Variance", las = 1, font = 1, cex = 1.3)
    }


    for(i in 2){
        par(mar = c(5,0,2,13))
        plot_marginal(marg_list_M1[[i]], expression(PV[max]), ylim = c(0,7), xexpan = c(-1.5, 1.5), cex.lab = 2, cex.axis = 1.5, plotname = FALSE)
        text(marg_list_M1[[i]]$Freq, x = 2, y = marg_list_M1[[i]]$obs, cex = 2)
        text("N", x = 2, y = max(marg_list_M1[[i]]$obs)+1, font = 2, cex = 2)
        mtext(at = 7, side = 2, "B)", las = 1, font = 2.5, cex = 2.5, adj = 0.7, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Decreased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  expression("variance along"~bold(g[max])), las = 1, font = 3, cex = 1.3)
        text(y = 7, x = 1, expression(bold("Increased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, expression("variance along"~bold(g[max])), las = 1, font = 3, cex = 1.3)
    }

    #dev.off()

Figure 3

    #pdf(file = "Fig3.pdf", height = 7.055555, width = 17.513889)
    par(mfrow = c(1,2))
    marg_list_M2 <- list(angle_within_marg_means, marg_evolve_max, marg_evolve_angle)
    marg_list_M2 <- lapply(marg_list_M2, function(x) cbind(x, name = c("Non-novel", "Novel", "Full-sib", "Half-sib")))
    marg_list_M2 <- lapply(marg_list_M2, function(x) cbind(x, obs = c(4,5,1,2)))


    labels = c("Breeding Design", "Environment Type")

    for(i in 2){
        par(mar = c(5,13,2,0))
        plot_marginal(marg_list_M2[[i]], expression(pi[e]), ylim = c(0,7), xexpan = c(-2.5, 1.5), cex.lab = 2, cex.axis = 1.5)
        #text(marg_list_M2[[i]]$Freq, x = 2, y = marg_list_M2[[i]]$obs, cex = 2)
        mtext(at = 3, side = 2, labels[1], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 6, side = 2, labels[2], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 7, side = 2, "A)", las = 1, font = 2.5, cex = 2.5, adj = 0.9, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Low")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  "Genetic Variance", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1, expression(bold("High")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, "Genetic Variance", las = 1, font = 1, cex = 1.3)
    }

    for(i in 3){
        par(mar = c(5,0,2,13))
        plot_marginal(marg_list_M2[[i]], expression(theta[e]), ylim = c(0,7), xexpan = c(-1.5, 2.5), cex.lab = 2, cex.axis = 1.5, plotname = FALSE)
        text(marg_list_M2[[i]]$Freq, x = 2, y = marg_list_M2[[i]]$obs, cex = 2)
        text("N", x = 2, y = max(marg_list_M2[[i]]$obs)+1, font = 2, cex = 2)
        mtext(at = 7, side = 2, "B)", las = 1, font = 2.5, cex = 2.5, adj = 0.7, padj = -0.8)
        text(y = 7, x = -1.2, expression(bold("G More Aligned")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1.2, "with plasticity vector", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1.1, expression(bold("G Less Aligned")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1.1, "with plasticity vector", las = 1, font = 1, cex = 1.3)
    }

    #dev.off()

Figure 4 Code

Note that this is not the complete figure. It was modified in Adobe Illustrator.

            Gn_e <- subset(evolve_max_data, environment == "Novel")
            Ga_e <- subset(evolve_max_data, environment == "Non-novel")

            Gn_angle_e <- subset(evolve_angle_data, environment == "Novel")
            Ga_angle_e <- subset(evolve_angle_data, environment == "Non-novel")

        par(mfrow = c(1,2), mar = c(6,6,1,1))
        adj = -0.15
        padj = 0.70
        #emax
        plot(Gn_e$ES ~ Ga_e$ES, ylab = expression(bolditalic(pi[e])~" Novel"), xlab = expression(italic(pi[e])~" Non-Novel"), pch = 16, cex = 1, ylim = c(-5, 16), xlim = c(-5, 16), cex.lab = 1.5, las = 1)
        abline(a = 0, b = 1, lty = 2, lwd = 1.5)
        abline(h = 0, col = "gray", lty = 2, lwd = 1.5)
        abline(v = 0, col = "gray", lty = 2, lwd = 1.5)
        
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = Ga_e$ES, y1 = (Gn_e$ES+(1.96*sqrt(Gn_e$S_var))), length = 0, angle = 90)
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = Ga_e$ES, y1 = (Gn_e$ES-(1.96*sqrt(Gn_e$S_var))), length = 0, angle = 90)

        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = (Ga_e$ES+(1.96*sqrt(Ga_e$S_var))), y1 = Gn_e$ES, length = 0, angle = 90)
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = (Ga_e$ES-(1.96*sqrt(Ga_e$S_var))), y1 = Gn_e$ES, length = 0, angle = 90)
        #text(x = Ga_e$ES+0.2, y = Gn_e$ES+0.2, Gn_e$study, cex = 0.2)
        mtext(expression(bolditalic("A)")), adj = adj, padj = padj,cex = 2)


        adj = -0.18
        padj = 0.70
        # e angle
        plot(Gn_angle_e$ES ~ Ga_angle_e$ES, ylab = expression(bolditalic(theta[e])~" Novel"), xlab = expression(italic(theta[e])~" Non-Novel"), pch = 16, cex = 1, ylim = c(-8.5, 6), xlim = c(-8.5, 6), cex.lab = 1.5, las = 1)
        abline(a = 0, b = 1, lty = 2, lwd = 1.5)
        abline(h = 0, col = "gray", lty = 2, lwd = 1.5)
        abline(v = 0, col = "gray", lty = 2, lwd = 1.5)

        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = Ga_angle_e$ES, y1 = (Gn_angle_e$ES+(1.96*sqrt(Gn_angle_e$S_var))), length = 0, angle = 90)
        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = Ga_angle_e$ES, y1 = (Gn_angle_e$ES-(1.96*sqrt(Gn_angle_e$S_var))), length = 0, angle = 90)

        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = (Ga_angle_e$ES+(1.96*sqrt(Ga_angle_e$S_var))), y1 = Gn_angle_e$ES, length = 0, angle = 90)
        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = (Ga_angle_e$ES-(1.96*sqrt(Ga_angle_e$S_var))), y1 = Gn_angle_e$ES, length = 0, angle = 90)

        #text(x = Ga_angle_e$ES+0.2, y = Gn_angle_e$ES+0.2, Gn_angle_e$study, cex = 0.2)
        mtext(expression(bolditalic("B)")), adj = adj, padj = padj, cex = 2)

        # Identify values where Gn_e CI's overlap the CI of the Ga_e


    dev.copy2pdf(file = "figure3_standardised10Mar.pdf", useDingbats=FALSE)

Figure 4 – A) The total amount of genetic variation in the direction of the plastic response vector as a proportion of the maximum amount of variation in any direction, π_e. Smaller values on the axes (i.e., negative numbers on axis) correspond to situations where genetic variation along the plasticity is low (‘Low’ label) whereas larger positive values correspond to situations where genetic variation along the plasticity vector is high (‘High’ label) Insets show schematic multivariate ellipses for the Non-novel (‘green’) and novel (‘orange’). Dashed arrows are the plastic response vector (the difference between centroids of multi-variate phenotype space across environments). Arrow length corresponds to ‘smaller’ (negative numbers on axis) or ‘larger’ (positive numbers on axis) amounts of genetic variation along the plasticity vector for the Non-novel (‘green’ text) and novel (‘orange’ text) environments. B) the The angle between the plastic response vector and g_max (θ_e) in novel and Non-novel environments across all studies. Labels correspond to values on the axes where the plasticity vector does not align with the major axis of genetic variation, g_max , (‘Not Aligned’) or where g_max aligns with the plasticity vector (‘Aligned’) in both the Non-novel (‘green’ text) and novel (‘orange’ text) environments. Insets show schematic angles between the plastic response vector (‘dashed arrow’) and g_max in the Non-novel (‘green’) and novel (‘orange’) environments. Axes are on transformed scale to facilitate comparisons with analyses. Large negative values (e.g., -8) indicate that the plastic response vector and g_max is fully aligned (i.e., ~0 degrees), whereas large positive values (e.g., 8) indicate that the plastic response vector and g_maxis not aligned (i.e., ~90). Values of 0 indicate that the plastic response vector and g_max are 45 degrees to one another. For both A) and B), studies that deviate from the expected 1:1 (dashed line) under the assumption of a constant G are colored, and indicate changes in the alignment between G and the plastic response vector across environments. ‘Green’ effects are those where the plasticity vector and G aligns better in the Non-novel environment relative to the novel environment. In contrast, ‘orange’ effects are effects where the plasticity vector and G aligns better with the novel environment relative to the Non-novel environment.

Figure 5 Code

    pdf(height = 9, width = 8)
    titles <- c("A) Kasule 1991", "B) Blanckenhorn & Heyland 2005", "C) Via & Connor 1995","D) Kause & Morin 2001")  

    # locations of the data in the dscrp and mats lists.
    findmeans <- c(6,8,28,20)
    findmats <- rbind(c(22,21,24,23),c(29,30,31,32),c(110,109,112,111),c(78,77,80,79))

    #par(mfrow=c(2,2))
    for each of the for studies:
    for(t in 1:4){
     
      # extract matrices
      mat <- list()
      mat[[1]] <- mats[[findmats[t,1]]]
      mat[[2]] <- mats[[findmats[t,2]]]
      mat[[3]] <- mats[[findmats[t,3]]]
      mat[[4]] <- mats[[findmats[t,4]]]
      
      # extract mean values
      mns <- array(NA,dim=c(2,4))
      mns[,1] <- dscrp[[findmeans[t]]]$mean[1:2]
      mns[,2] <- dscrp[[findmeans[t]]]$mean[1:2]
      mns[,3] <- dscrp[[findmeans[t]]]$mean[3:4]
      mns[,4] <- dscrp[[findmeans[t]]]$mean[3:4]#
      # colors for each matrix
      cols <- c("#348B08A0","#54AB28A0","#DE8601A0","#FEA621A0")#
      # calculate eigenvalues and vectors
      eig <- lapply(mat,eigen)
      
      # calculate coordinates of vectors
      crs <- lapply(eig,function(x)x$vectors*matrix(sqrt(x$values),nrow=2,ncol=2,byrow =TRUE))
      
      # calculate the plotting range
      xlim <- ylim <- list()
      for(i in 1:4){
        x1 <- mns[,c(i,i)]+1.5*crs[[i]]
        x2 <- mns[,c(i,i)]-1.5*crs[[i]]
        xlim[[i]] <- c(x1[1,],x2[1,])
        ylim[[i]] <- c(x1[2,],x2[2,])  
        }
      xlim <- range(unlist(xlim))
      ylim <- range(unlist(ylim))#
      # draw plotting plane 
      plot(NA,xlim=xlim,ylim=ylim,xlab="trait 1",ylab="trait 2", main=titles[t])#
      # add ellipses and their axes for the G and P-matrices
      for(i in 1:4){
        ell <- car::ellipse(mns[,i], as.matrix(mat[[i]]), draw=FALSE, radius=1, center.pch=16, center.cex=1)
        polygon(ell[,1],ell[,2], col=cols[i])
        for(j in 1:2){
          segments(mns[1,i],mns[2,i],mns[1,i]+crs[[i]][1,j],mns[2,i]+crs[[i]][2,j])
          segments(mns[1,i],mns[2,i],mns[1,i]-crs[[i]][1,j],mns[2,i]-crs[[i]][2,j])    
          }
      }
    ## add plasticity vector
    segments(mns[1,1],mns[2,1],mns[1,3],mns[2,3],lwd=3, lty=2) 
    Arrows(mns[1,1],mns[2,1],mns[1,3],mns[2,3],lwd=1, segment ="false" ,arr.type="triangle", arr.width=0.3 , arr.length=0.4, arr.adj = 1)   
    }
    dev.off()

Figure 5 – Examples of studies (with two trait matrices) with high (A, B) or low (C, D) alignment between the direction of plastic responses and standing genetic variation. Arrows are plastic response vectors; large ellipses are P matrices and smaller ellipses are G matrices. Colors correspond to the Non-novel (i.e., ‘green’) and novel (i.e., ‘orange’) environments. A) Kasule (25) used the cotton stainer bug (Dysdercus fasciatus) to estimate quantitative genetic variation in fecundity and adult size under stressful conditions using a half-sib breeding design. B) Blanckenhorn & Heyland (26) used a full-sib breeding design to estimate quantitative genetic variation for hind limb size and development time in the yellow dung fly (Scathophaga stercoraria) under stressful conditions. C) Via & Conner (27) studied the flour beetle (Tribolium castaneum) using a half-sib breeding design to estimate quantitative genetic variation for pupal weight and development time under novel (not stressful) conditions. D) Kause & Morin (28) studied the sawfly (Priophorus pallipes), and estimated quantitative genetic parameters for development time and body size using a full-sib breeding design (broad-sense estimates) under stressful environmental conditions.

6.0 – Additional Supplemental Tables

6.1 – Table S2

Table S2 – Model estimates for (1) standardized differences in total trait genetic variance (SDV) between Non-novel and novel environments, (2) the change in proportion of variance along g_max (PV_max) between Non-novel and novel environments, and (3) the angle between g_max (θ_G) between Non-novel and novel environments. 95% confidence intervals are provided. Intercept values are the average effect size for a study with an average number of traits conducting a full-sib breeding design in a benign environment. All other parameters are contrasts from this mean. Estimates in boldface indicate significant effects (P < 0.05).

	Est. (1)	L 95 %	U 95 %	Est. (2)	L 95 %	U 95 %	Est. (3)	L 95 %	U 95 %
Intercept	0.623	0.086	1.159	0.465	-0.501	1.431	-1.97	-3.015	-0.925
Half-sib Breeding Design	-1.092	-1.749	-0.434	-1.475	-2.676	-0.275	-1.245	-2.551	0.062
Stressful Environment	-0.012	-0.634	0.610	0.751	-0.387	1.889	1.51	0.269	2.752
Number of Traits	-0.090	-0.358	0.179	0.454	-0.044	0.952	0.89	0.353	1.427

6.2 – Table S3

Table S3 – Estimates of the (1) the angle between g_max and p_max within environments, (2) total amount of genetic variance in the direction of plastic response vector (as a proportion of g_max), π_e, and (3) the angle between the plastic response vector and g_max, θ_e, in in both the novel and Non-novel environments. 95% confidence intervals are provided. Intercept values are the average effect size for a study with an average number of traits conducting a full-sib breeding design in an Non-novel environment. All other parameters are contrasts from this mean. Estimates in boldface indicate significant effects (P < 0.05).

	Est. (1)	L 95 %	U 95 %	Est. (2)	L 95 %	U 95 %	Est. (3)	L 95 %	U 95 %
Intercept	-2.083	-3.781	-0.385	1.274	-0.353	2.900	-0.525	-1.569	0.520
Novel Environment	0.094	-0.251	0.439	0.088	-0.420	0.597	-0.102	-0.458	0.254
Half-sib Breeding Design	1.080	-0.995	3.155	-0.756	-2.689	1.177	0.449	-0.790	1.689
Number of Traits	0.080	-1.876	2.036	-1.074	-1.875	-0.273	0.733	0.213	1.254

7.0 – Functions Used Throughout

Here we provide the details of the various functions that we use throughout so that readers can have a look at what they are doing in more detail.

## Comparing G and P matrices Functions

# Write files from a list
    writeFiles <- function(listFiles, dir="./data/", add = ".csv"){
        names <- paste0(dir, names(listFiles), add)
        
        for(i in 1:length(names)){
            write.csv(listFiles[[i]], file = names[i])
        }
    }

# Function will go into the data folder and grab all the relevant matrix files ("matrix")or the descriptive data ("D"). It will return a list of folder names with the files within these contents, ignoring irrelevant folders. 

    files <- function(dir = "./data/", file_pattern = "[WRs][BR]", type = c("D", "matrix")){
        type <- match.arg(type)
         tmp <- paste0(dir, list.files(dir)[grep(file_pattern, list.files(dir))])

        if(type == "D"){
            file_names <- lapply(tmp, function(x) paste0(x, "/",list.files(x)[grep("Data", list.files(x))]))
        }

        if(type == "matrix"){
            file_names <- lapply(tmp, function(x) paste0(x, "/",list.files(x)[grep("n=", list.files(x))]))
        }
        names(file_names) <- tmp
        return(file_names)
    }


# Importing matrices
    
    import_mat <- function(list){
            tmp <- as.character(unlist(list))
            mat <- suppressWarnings(lapply(tmp, function(x) read.csv(file = x, header = TRUE)))
            names(mat) <- tmp

            return(mat)
    }


## Function for extracting the sample sizes from the file names. Argument takes the string name for the various matrices and extracts the sample sizes, returning a numeric vector.

    match_n <- function(string){
            m <- regexpr("=[0-9]*(_|.)", string)
            tmp <- regmatches(string, m)
            tmp2 <- gsub("=","", tmp)
            tmp3 <- as.numeric(gsub("_","", tmp2))
        return(tmp3)
    }

## Function extracts the study names of the folders from each of the matrices   
    extract_matrix_studynames <- function(string, pattern = "[WRs][BR][0-9][0-9][0-9].*/"){
            m <- regexpr(pattern, string)
            tmp <- regmatches(string, m)
            tmp2 <- gsub("/","", tmp)
            return(tmp2)
    }

# Extract Environments
    
    extract_matrix_env <- function(string, pattern = "/[ABNSU].", type = c("full", "reduced")){
            m <- regexpr(pattern, string)
            tmp <- regmatches(string, m)
            tmp2.1 <- gsub("/","", tmp)
            tmp2.2 <- gsub("_","", tmp2.1)
            tmp2.3 <- gsub("[0-9]","", tmp2.2)

            if(type == "full"){
                tmp3 <- gsub("_","", tmp2.1)
                return(tmp3)
            }

            if(type == "reduced"){
                return(tmp2.3)  
            }
    }


# Extract whether matrix is a G or P
    extract_matrix_type <- function(string, pattern = "_[GP].csv"){
                m <- regexpr(pattern, string)
                tmp <- regmatches(string, m)
                tmp2 <- gsub("_","", tmp)
                tmp3 <- gsub(".csv","", tmp2)
                return(tmp3)

    }


    ## Simulation function. This will use Monte Carlo simulations to generate various statistics from a given matrix, assuming the traits follow a multivariate normal distribution. A number of different statistics can be chosen, including the total variance in a matrix Vt or "totalVar", which can be used to then look at the difference between two matrices; the proportion of variation along g/p max "varAlong", or it can even extract the full set of eigen values from a matrix ("eigen") or only the G/P max of a matrix ("max"). This works by simply taking a matrix, with it's corresponding sample size and assuming the traits follow a multivariate normal distribution. TThe matrix is used to simulate 1000 matrices according to this MVN distribution, generating the various statistics from each simulated dataset that follows this matrix structure. From these we can get a fairly robust estimate of sampling variance for each matrix and this entire simulated distribution can be used to generate effect sizes.
# x = means
     sim <- function(matrix, x, n, sims = 1000, cor = FALSE, type = c("totalVar", "varAlong", "raw_matrix"), logit = TRUE, matrix_fix = c("shrinkage", "bend"), mu = rep(0, nrow(matrix))){
              
                if(corpcor::is.positive.definite(matrix) == FALSE){

                    if(matrix_fix == "shrinkage"){
                        matrix <- corpcor::cov.shrink(matrix)
                    }
                    
                    if(matrix_fix == "bend"){
                        matrix <- make.positive.definite(matrix)
                    }       
                }

                    eigens <- list()
                    
                    if(type == "totalVar"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          eigens[[i]] <- sum(eigen(cov)$values)
                        }
                    }

                    if(type == "varAlong"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          
                          if(logit){
                            p <- (eigen(cov)$values[1]) / sum(eigen(cov)$values)
                            eigens[[i]] <- log(p / (1-p))       
                          } else{
                            eigens[[i]] <- (eigen(cov)$values[1]) / sum(eigen(cov)$values)      
                          }
                        }
                    }

                    if(type == "raw_matrix"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          eigens[[i]] <- cov
                        }
                    }

                    if(type == "raw_matrix"){
                        return(eigens)
                    }else{
                        return(do.call(rbind,eigens))
                    }
     }


# Function for matrix bending by Reinder

    make.positive.definite <- function(x, tol=1e-6){
      eig <- eigen(x, symmetric=TRUE)
      rtol <- tol * eig$values[1]
      if(min(eig$values) < rtol) {
        vals <- eig$values
        vals[vals < rtol] <- rtol
        srev <- eig$vectors %*% (vals * t(eig$vectors))
        dimnames(srev) <- dimnames(x)
        return(srev)
      } else {
        return(x)
      }
    }

# Angle between two vectors a and b
    angle <- function(a,b){
      
      angle <- acos( sum(a*b) / ( sqrt(sum(a * a)) * sqrt(sum(b * b)) ) ) * (180/pi)

      angtmp <- ifelse(angle > 90, 180-angle, angle)

      # Logit transformation of the angle to normalise the distribution.
      ang <- log( (angtmp/90) / (1 - (angtmp/90) ) )
      
      return(ang)

    }

## Use function below to calculate differences between effects within a study and one environment between G and P. Note that P is subtracted from G. NOTE::: "angle" can ONLY be used with a simulated effect of "max". p = the proportion of vectors, rounded to compare for angles and similarity. x = the columns from the metadata frame. Stats is the simulation object.
    
#type = "SMD"
#p = 0.5
#environment = "across"
#stats = sim_object
#x = metadata
    
    compare_within_env <- function(x, stats, type = c("raw_diff", "SMD", "angle"), environment = c("within", "across")){
            
            cols <- as.vector(x$col)   

        if(type == "raw_diff"){
            col_compute <- data.frame( SMD = (stats[, cols][,2] - stats[, cols][,1]) )
        }

        if(type == "SMD"){
            col_compute <- data.frame( SMD = (stats[, cols][,2] - stats[, cols][,1]) / ( (stats[, cols][,1] + stats[, cols][,2] ) / 2) )
        }


        if(type == "angle"){
                # NOTE::: this can ONLY be used with simulated matrices
                matrix1 <- stats[, cols][,1]
                matrix2 <- stats[, cols][,2]
        
                eig1 <- lapply(matrix1, function(x) eigen(x)$vectors[,1])
                eig2 <- lapply(matrix2, function(x) eigen(x)$vectors[,1])

                # Compute angle
                angle <- mapply(function(a,b) angle(a,b), a = eig1, b = eig2)   
                
                col_compute <- as.data.frame(angle)
            }

        if(environment == "within"){
            colnames(col_compute) <- unique(x$nam.env)
        } else{
            colnames(col_compute) <- unique(x$nam.type)
        }
        
            return(col_compute)
        
    }


# After simulating and creating effects sizes, create the data frame. Change the type of effect size that you wish to generate using the type argument. The choice of this really depends on the specific simulated product that you get. The choices are: c("raw_diff", "SMD", "angle_krzanowski", "angle", "similarity_krzanowski"). raw_diff is simply the subtraction of TOTAL variance of the matrix and the SMD is the standardised mean difference, where the difference in variance is standardised by the combined total variance of each of the matrices. This is important to make sure that G and P are put on comparative scales for between environment comparisons of each matrix.

create_data <- function(sim_object, type = c("SMD", "angle"), environment = c("within", "across"), plot = c("yes", "no")){
    
              nam <- extract_matrix_studynames(colnames(sim_object))
              env <- extract_matrix_env(colnames(sim_object), type = "reduced")
        type_name <- extract_matrix_type(colnames(sim_object))

    # Create a metadata file to summarize and check that everything so far has gone smoothly and as expected.
        
        metadata <- data.frame( file = colnames(sim_object), 
                                 nam = nam, 
                                 env = env, 
                                type = type_name, 
                                   n = n, 
                                 col = 1:length(colnames(sim_object)), 
                                 stringsAsFactors = FALSE )

            metadata$nam.env  <- paste0(metadata$nam, ".", metadata$env)
            metadata$nam.type <- paste0(metadata$nam, ".", metadata$type)

    ## Calculate the difference between G & P within an environment
    
    ## Make sure that every study is arranged by environment and G/P in the same way so that calculations between columns are done correctly. Or, if interested in the GP comparison across environments then we can arrange slightly differently

    # Split metadata on study and environment
        if(environment == "within"){

            split_within <-  split(metadata, metadata$nam.env, drop = TRUE)
            split_within <-  lapply(split_within, function(x) dplyr::arrange(x, type))

        } else {

            split_within <- split(metadata, metadata$nam.type, drop = TRUE)
            split_within <-  lapply(split_within, function(x) dplyr::arrange(x, env))

        }
        

    # Calculate for all environments 
        
        effect_size <- do.call(cbind, lapply(split_within, function(x) compare_within_env(x, sim_object, type = type, environment = environment)))


    # From the simulated distributions, calculate the mean and the sampling variance (i.e., the variance of the distribution)
           
              mean <- sapply(effect_size, mean)
                Vi <- sapply(effect_size, var)
    
    #Create a dataframe
    
    if(environment == "within"){
        dat   <- data.table(study = gsub("[.][SNAB]", "", names(mean)), 
                            env = gsub(".+[.]", "", names(mean)), 
                            ES = mean, 
                            S_var = Vi )
    } else{
        dat   <- data.table(study = gsub("[.][GP]", "", names(mean)), 
                            type = gsub(".+[.]", "", names(mean)), 
                            ES = mean, 
                            S_var = Vi )
    }
    
    # Final data for differences in g/pmax within an environment
        
    if(environment == "within"){
        data <- dplyr::arrange(dat, study, env)
    } else{
        data <- dplyr::arrange(dat, study, type)
    }

    if(plot == "yes"){
            tmp <- effect_size %>% gather()
            fig <- ggplot(tmp, aes(value)) + 
            geom_histogram() + 
            facet_wrap(~key, scales = 'free_x') + theme(axis.text.x=element_text(size=rel(0.5))) + labs(x = "Effect Size")
            return(fig)
        } else{
            return(data)    
        }

}
    


## Evolvability Specific Functions

evolveability_along_plastic_vector <- function(G_sim, cen.diff, compare = "max"){
        
                colnames(G_sim) <- extract_matrix_studynames(colnames(G_sim))
            names_cen.diff_evolve <- names(cen.diff)

            # check that names match
            if(sum(colnames(G_sim)!=names_cen.diff_evolve)>0) {stop("Sorry, you have different matrices and response vectors") }
            
            evol <- c()
            
            for(j in names_cen.diff_evolve){
                        
                        G <- G_sim[,j]
                    r_vec <- cen.diff[[j]]
                
                    tmp <- do.call("rbind", lapply(G, function(x) project.matrix(x, r_vec, compare = compare)))
                    colnames(tmp) <- j
                    
                    evol <- cbind(evol, tmp)
            }

            return(as.data.frame(evol))
}


evolveability_angle_along_plastic_vector <- function(G_sim, cen.diff){
        
                colnames(G_sim) <- extract_matrix_studynames(colnames(G_sim))
            names_cen.diff_evolve <- names(cen.diff)

            # check that names match
            if(sum(colnames(G_sim)!=names_cen.diff_evolve)>0) {stop("Sorry, you have different matrices and response vectors") }
            
            evol <- c()
            
            for(j in names_cen.diff_evolve){
                
                    # Grab the same study for response vector and G     
                        G <- G_sim[,j]
                    r_vec <- cen.diff[[j]]
                    
                    # Calculate the angle between gmax and response vector
                    tmp <- do.call("rbind", lapply(G, function(x) angle(eigen(x)$vectors[,1], r_vec)))
                    colnames(tmp) <- j
                    
                    evol <- cbind(evol, tmp)
            }

            return(as.data.frame(evol))
}

create_data_evolve <- function(sim_object, meta_evolve, env = c("N", "S")){
    
        columns <- c("authors", "year", "study", "env", "type", "n", "trait_num", "genus", "species", "design", "group", "name.type", "name.env")
    trim_meta <- meta_evolve[meta_evolve$type == "G", columns] 
    trim_meta <- trim_meta[trim_meta$env %in% env,]

            mean <- sapply(sim_object, mean)
             Vi <- sapply(sim_object, var)
             study <- names(mean)

             dat   <- data.table(study=study,  
                            ES = mean, 
                            S_var = Vi )

             data_fin<- left_join(trim_meta, dat, by = "study")

             return(data_fin)
}


combine_meta_data_evolve <- function(es_data, meta_sub){
            
                es_data <- unique(left_join(es_data, meta_sub, by = c("file")))

            # Fix up a couple columns
            es_data$species_full <- paste0(es_data$genus, "_", es_data$species)
            es_data$stdy <- substr(es_data$study, 1, 5)
            rownames(es_data) <- 1:nrow(es_data)
        return(es_data)
}


sim_G <- function(Gmatrix, n, sims = 1000, cor = FALSE, matrix_fix = c("shrinkage", "bend")){
                    
                    mu <- rep(0, nrow(Gmatrix))
                
                if(corpcor::is.positive.definite(Gmatrix) == FALSE){

                    if(matrix_fix == "shrinkage"){
                        Gmatrix <- corpcor::cov.shrink(Gmatrix)
                    }
                    
                    if(matrix_fix == "bend"){
                        Gmatrix <- make.positive.definite(Gmatrix)
                    }       
                }

                matrices <- list()
                    
                        for(i in 1:sims){
                          
                          d <- MASS::mvrnorm(n, mu, Gmatrix)
                          
                          if(cor == FALSE){
                            cov <- cov(d)
                          } else {
                            cov <- cor(d)
                          }
                          
                          matrices[[i]] <- cov
                        }       

                return(matrices)                    
}

centroid.diff <- function(a,b){
  # a: vector a
  # b: vector b
  dz <- b-a
  val <- sqrt(sum(dz^2))
  vec <- dz/val
  return(list(value=val,vector=vec))
}

# Function to project a covariance matrix onto a vector.

project.matrix <- function(mat,vec,compare=NA){
  # mat: covariance matrix that needs to be projected
  # vec: vector on which the covariance matrix needs to be projected
  # compare: to what does the amount of variance of mat on vec needs to be compared to:
  #   NA        : no comparison, output is the amount of variance on vec
  #   "max"     : as proportion of the first eigenvector of mat
  #   "maxcorr" : as proportion of the first eigenvector of mat, corrected for the expected value
  #   "mean"    : as proportional difference to mean value of all eigenvectors of mat
  #   "total"   : as proportion of the sum of the values of all eigenvectors of mat
  
  dims       <- length(vec)
  beta_prime <- matrix(vec,ncol=dims,nrow=dims)
  beta       <- matrix(vec,ncol=dims,nrow=dims,byrow=TRUE)
  abs_beta_2 <- sum(vec^2)
  pv <- as.numeric((t(vec)%*%mat%*%vec)/abs_beta_2)
  ev <- eigen(mat)$values
  if(is.na(compare)){
    return(pv)
  }else if(compare=="max"){
    return(log( (pv/ev[1]) / (1 - (pv/ev[1]))))
  }else if(compare=="maxcorr"){
    return(log( (pv/ev[1]) / (1 - (pv/ev[1])))-log( (mean(ev)/ev[1]) / (1 - (mean(ev)/ev[1]))))
  }else if(compare=="total"){
    return(pv/sum(ev))    
  }else{
    return((pv-mean(ev))/mean(ev))    
  }
}


## Adding some new functions for calculating marginalised means. 

# Relevels the data based on a factor in the model and a reference level of the factor
    re_level_dat <- function(data, factor, ref){
        data[, factor] <-   relevel(data[, factor], ref = ref)
        return(data)
    }

    #@example test <- re_level_dat(evolve_max_data, "design", ref = "half-sib")


## Extract coefficients from models
    extract_coefs <- function(model){
        mean <- model$b[1]
          se <- model$se[1]

          dat <- data.frame(mean = mean, se = se)

     return(dat)
    }

# Calculate marginal means for a single factor
    marg_mean <- function(coefs){
        mean_mar <- ( (coefs[1,"mean"]*coefs[1,"n"]) + (coefs[2,"mean"]*coefs[2,"n"]) ) / sum(coefs$n)
        return(mean_mar)
    }

    #@examples marg_mean(coefs)

#Calculate marginal SE's for a single factor

    marg_se <- function(coefs){
        n1 <- (coefs[1,"n"] / sum(coefs$n))
        n2 <- (coefs[2,"n"] / sum(coefs$n))

        mar_se  <- sqrt((coefs[1,"se"]^2)*(n1^2) + (coefs[2,"se"]^2)*(n2^2))
        return(mar_se)
    }

    #@example marg_se(coefs)

# Function for producing marginalised mean and se across a single factor level. To get estimates for a 2x2 factor design simply refit the model re-leveling the second factor (i.e., not the one used in the "marginal" function)
    
    marginal <- function(model, factor, data){
            
        # Get relevant data
                 ref <- levels(data[ , factor])
            n_levels <- table(data[ , factor])

        # Get first set of estimates from model
               coef1 <- cbind(extract_coefs(model), name = ref[1])

        # Update data, to get the new estimates
             dat_new <- re_level_dat(data, factor, ref = ref[2])

        # Refit model
            model_new <- update(model, data = dat_new)

        # Extract updated coefs
            coef2 <- cbind(extract_coefs(model_new), name = ref[2])

        # Concatenate the two rows
            coefs <- cbind(rbind(coef1, coef2), n = as.vector(n_levels))

        #Calculate the marginal mean and SE across the levels of "factor"
            marginal <- data.frame(mean = marg_mean(coefs),
                                   se   = marg_se(coefs)  )
        
        return(marginal)
    }

 #  @example{
#           full_A <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)#

#           marginal(full_A, "design", evolve_max_data) # Marginal mean of Non-novel across design#

#           full_N <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)#

#           marginal(full_N, "design", evolve_max_data) # Marginal mean of Novel across design
#   }


    plot_marginal <- function(table, effectsize, xexpan, plotname = TRUE, ...){
        plot(obs ~ mean, type = "n", las = 1, ylab = "", xlab = effectsize, data = table, bty = "n", yaxt = "n",  xlim = range(table$mean)+xexpan, ...)
        abline(v=0, lty = 2, lwd = 1.5)
        points(obs ~ mean, data = table, pch = 16, cex = 1.8)

        arrows(x0=table$mean, y0=table$obs, x1=table$mean + 1.96*(table$se), y1=table$obs, length = 0, angle = 90)
        arrows(x0=table$mean, y0=table$obs, x1=table$mean - 1.96*(table$se), y1=table$obs, length = 0, angle = 90)  

        if(plotname){
            mtext(table$name, at = table$obs, side = 2, cex = 1.3, las = 1)
        }
    }